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INTRODUCTION 


This  report  presents  the  results  of  an  investigation  for  representing 
the  transient  electromagnetic  coupling  through  a rectangular  aperture  by 
means  of  the  singularity  expansion  method. 

The  singularity  expansion  method,  which  was  introduced  by  Baum  in  1971 
(ref.  1),  has  been  subsequently  applied  to  many  scatterer  geometries. 

The  essence  of  thj  singularity  expansion  method  is  the  representing  of 
the  temporal  response  of  a body  in  terms  of  the  complex  natural  frequencies 
for  the  body. 

Taylor  et  al.  point  out  that  the  singularity  expansion  for  an  aperture 
in  an  infinite  perfectly  conducting  screen  can  be  determined  in  terms  of  that 
for  the  complementary  perfectly  conducting  plate  by  way  of  Babinet's 
principle  (ref.  2).  This  approach  was  taken  in  the  work  reported  here.  The 
remaining  discussion  is  directed  to  the  equivalent  problem  of  determining 
the  current  distributions  on  the  complementary  plate  geometry. 

I 

Rahraat  -Samii  and  Mittra  have  derived  a coupled  pair  of  Hallen-type 
integral  equations  governing  the  current  behavior  on  the  rectangular  plate 
(ref.  3).  The  work  reported  here  builds  on  their  work  by  generalizing  the 
integral  equations  and  solution  method  to  the  complex  frequency  plane  for  the 


1.  Baum,  C.  E.,  "On  the  Singularity  Expansion  Method  for  the  Solution  of 
Electromagnetic  Interaction  Problems,"  Interaction  Note  88,  Air  Force 
Weapons  Laboratory,  Kirtland  AFB,  NM,  December  1971. 

2.  Taylor,  C.  D.,  Crow,  T.  T.,  and  Chen,  K-T,  "On  the  Singularity  Expansion 

Method  Applied  co  Aperture  Penetration:  Part  I Theory,"  Interaction  Note 

134,  Air  Force  Weapons  Laboratory,  Kirtland  AFB,  NM,  May  1973. 

3.  Rahmat-Samii,  Y.  and  Mittra,  R.,  "Integral  Equation  Solution  and  RCS 
Computation  of  a Thin  Rectangular  Plate,"  Interaction  Note  156,  Air  Force 
Weapons  Laboratory,  Kirtland  AFB,  NM,  December  1973. 
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SEM  application.  The  same  method-o£-moments  formulation,  as  described  in 
(ref.  3),  is  used,  i.e.,  two-dimensional  pulse  expansion  functions  with 
collocation  testing. 

In  order  that  the  computation  time  be  practical  in  a problem  of  this 
complexity,  a great  deal  of  care  was  given  to  algorithmic  streamlining  in  the 
conduct  of  this  work.  The  streamlining  includes  maximum  exploitation  of 
geometric  symmetry,  organization  of  calculations  to  make  use  of  redun- 
dant terms  and  partial  terms  occurring  in  the  calculation,  and  direct 
algorithmic  exploitation  of  matrix  sparseneaa.  The  end  result  is  a highly 
efficient  computer  code.  Key  features  of  the  algorithms  are  dis- 
cussed in  this  report. 

The  pules  expansion  appears  to  be  inadequate  in  accurately  modeling  the 
rectangular  plate.  The  difficulty,  which  relates  to  representing  the  actual 
size  of  the  plate,  is  demonstrated  and  discussed  herein.  Remedies  for  the 
problem  are  suggested,  but  they  are  outside  the  scope  of  the  present  work. 

Ry  holding  the  zoning  scheme  invariant  while  the  aspect  ratio  of  the 
plate  was  changed,  self-consistent  pole  trajectories  for  tha  four  fundamental 
modes  were  determined.  For  the  reasons  cited  above,  these  poles  cannot  claim 
to  be  the  exact  poles  for  the  body.  They  are,  however,  indicative  of  the 
trends  in  the  pole  behavior  for  the  plate  under  change  in  aspect  ratio.  These 
results  are  reported  and  discussed  in  this  context. 
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SECTION  II 

THIN-PLATE  INTEGRAL  EQUATION  FORMULATION 
FOR  COMPLEX  WAVENUMBER 


Rahmat-Samii  and  Mittra  (ref.  3)  give  an  integral  equation  formulation 

for  the  rectangular  plate  subject  to  time-harmonic  excitation.  Their  results 

may  be  directly  extended  to  the  complex  wavenumber  case.  That  is,  for  the 

geometry  in  Figure  1 with  exp[st]  time  dependence,  s = a + jw  complex,  and  an 

incident  plane-wave  magnetic  field  component 

H*=[H^u  + H * u + H * u ] exp[j(k  x + k y + k z)]  the  following 

ox  x oy  y oz  z r J x yJz 

coupled  integral  equations  result: 


L/2  v/2 


-L/2  -w/2 


rfJx  (x*y) 

/w 


K(x,y|x',y')  dx'  dy*  = ® lexp[j(kx  x + k y)  ] 


z 1 -Ho 


+ 1 J [ exp[j(n  + 1)<(^  Jn+l^*sp/c) 


The  kernel  is  given  by 


kr  exp[j(n  " 1H]  Jn-1  (-s?/c>l 


K(x,y  | x',y')  = exp[-sR/c]/R 


R = [(x  - x')2  + (y  - y')2]J/2 

The  J (x,y)  and  J (x,jr)  denote  the  respective  x and  y components  of  current 
x y 

on  the  plate;  J (c)  denotes  the  Bessel  function  of  the  first  kind;  the  C are 
n n 

unknown  constants;  c is  the  velocity  of  light;  and  (p,4>)  are  the  polar  coor- 
dinates for  the  point  (x,y)  on  the  plate.  Equation  (1)  holds  for 
xe  (-L/2,L/2)  and  ye  (-w/2, w/2),  and  z = 0. 


It  is  pointed  out  that  the  two  integral  equations  represented  by  (1)  are 


coupled  through  the  in  the  summation  in  the  right-hand  side.  This  sum- 
* 

matlon  is  simply  a Bessel  function  expansion  of  the  homogeneous  solution  to 
the  wave  equation  which  occurs  in  the  derivation  of  (1) . Details  of  arriving 
at  this  expansion  are  found  in  (ref.  3).  The  pair  of  integral  equations  (1) 
is  complete  in  the  sense  that  no  additional  constraints  are  needed  to  cor- 
rectly specify  the  currents.  It  is  noteworthy,  however,  that  current 
solutions  to  (1)  satisfy  the  Meixner’s  edge  condition  (ref.  4).  Namely, 

T Ti/I  /■)  _ ..1  J 1/2  ^ 


Jx[±(L/2  - d),  y]  -►  d 

Jy[±(L/2  - d),  y]  d’ 

Jxtx,±(w/2  - d)]  d' 

Jy[x,±(w/2  - d)]  -►  d1 


d 0 


The  ability  of  a numerical  solution  to  approximate  the  behavior  of  eqn.  (3)  is 
a key  point  in  a subsequent  discussion. 


Bladel,  j.  Van,  Electromagnetic  Fields,  McGraw-Hill,  New  York,  pp.  385-387, 
1964. 
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SECTION  III 


SYMMETRY  CONDITIONS  FOR  THE  NATURAL  MODE  CURRENTS 


The  natural  frequencies  of  (1)  occur  when  the  complex  frequency  s is 

such  that  there  are  non-trivial  and  and  the  accompanying  Cr 

satisfying  (1)  for  H*  ■ 0.  Such  J and  J solutions  are  natural  mode 

x y 

current  solutions  for  the  rectangular  plate,  and  the  concomitant  value  of  s 
is  a pole  of  the  plate.  The  vanishing  of  incident  wave  dependence  gives 
rise  to  symmetry  in  the  integral  equations.  By  discerning  the  symmetry 
relations  a priori  and  bringing  them  to  bear  upon  solution  procedures,  one 
gains  significant  computational  savings  in  the  numerical  solution  for  poles 
and  natural  modes.  These  symmetry  relations  are  explored  in  this  section. 

The  excitation-free  form  of  (1)  is 
L/2  w/2  . r 

J J Jx  K » y I x * * y ' ) dx'  dy'  = l cn(jn+1  exP[j(n  + i)<j>]  Jn+1  (-sP/c) 

-L/2  -w/2 


.n-1 


+ j exp[  j (n— 1)  ] Jn_1  (-sp/c) 


•} 


(4a) 


and 
L/2  w/2 


C f °°  / 

J J Jy  K(x»ylx'»y'>  dx'  dy'  = "T  l cn(  j"+1  expt^n+D<j>]  Jn+1  c-sP/c) 

n-i 


-L/2  -w/2 


- jn  1 exp[j  (n  - 1) d>]  J, 


By  using  the  symmetry  of  the  Bessel  function  with  respect  to  order,  expanding 
the  exponentials  by  way  of  Euler's  identity,  and  appropriately  adjusting  the 
indices,  one  arrives  at  the  following  equation  after  some  manipulation. 
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L/2  w/2 

/ / Jx  K dlt'  “y' 

-L/2  -w/2 


ttc 

s 


and 


: I fjn+1  dn  [cos(n  + 1)4>  J^_1  (— sp / c)  - cos(n  - 1)4>  J . (-sp/c)  ] 
n-0  ’ 

- jnd“  [sin(n  + 1)4>  J^C-sp/c)  - sin(n  - 1)<|>  J^O-sp/c)^  ( 


L/2  w/2 

J J JyK  dx'  dy* 

-L/2  -w/2 


[sln(n  + 1)<})  J^C-sp/c)  + sin(n  - 1)4>  J^C-sp/c)] 


: 1 6n+1  < 

n-0  V. 

+ jndn  [cos(n+  1)4>  J^C-sp/c)  + un  l cos(n  - 1)4>  Jnl(-sp/c)ij 


(5b) 


where 


and 


u -A’  n ^ 0 

n \p,  n < 0 


It  is  noted  that  the  d*  multiply  terras  containing  cosine  functions  in  the  Jx 
equation,  while  they  multiply  terras  containing  sine  functions  in  the 
equation.  The  situation  is  reversed  for  the  d^. 

Because  of  the  symmetry  properties  of  the  kernel,  the  integral  operator 
on  the  left-hand  sides  of  (5)  produces  a function  whose  symmetry  character  is 
Identical  to  that  of  the  current  on  which  it  operates.  Then,  for  a given 

4* 

current  symmetry,  only  part  of  the  d^  on  the  right-hand  side  may  be  non- 
zero because  of  the  symmetries  possessed  by  the  trigonometric  terms.  Thus, 
the  respective  symmetries  for  and  J , which  are  compatible,  and  the 
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surviving  terms  in  the  right-side  series  may  be  discerned  by  1)  postulating  a 
symmetry  for  J^,  2)  determining  from  (5a)  which  right-hand  side  terms  survive  so 
as  to  be  compatible  with  the  symmetry,  3)  observing  in  (5b)  the  variation 
which  terms  have  non-zero  coefficients,  and  4)  determining  the  Jy  symmetry 
conditions  compatible  with  the  postulated  symmetry  conditions. 

For  example,  if  is  symmetric  with  respect  to  the  y axis  and  anti- 
symmetric with  respect  to  the  x axis,  orly  sin(n  + 1) 4>  terms  with  n even  are 
compatible  in  (5*»).  Thus,  cnly  dn  , n even,  may  be  non-zero.  In  the  right- 
hand  side  of  (5b),  the  coefficients  multiply  cos(n  + 1)<J>  terms  with  n even. 

These  cosines  sum  to  functions  which  are  antisymmetric  with  respect  to  the 
y axis  and  symmetric  with  respect  to  the  x axis.  Stated  mathematically,  if 


Jx(*,y)  “ Jx(-x»y) 

(6a) 

anti 

Jx(x,y)  ° -Jx(x,-y) 

(6b) 

then 

d+  « 0,  for  all  n, 

(6c) 

d“»  0,  n odd, 

n 

(6d) 

and 

Jy(x,y)  - -Jy(-x,y) 

(6e) 

Jy(x,y)  - Jy(x,-y) 

(6f ) 

These  vector  symmetries  are  in  accord  with  the  general  symmetry  relations 
given  by  Baum  (ref.  5).  The  information  in  (6)  may  be  used  to  reduce  the 
complexity  of  the  integral  equations  (4),  viz.,  by  (6a,b,e,f)  the  range  of 


each  integration  may  be  halved  while  by  (6c, d)  the  zero  terms  of  the  right- 
hand  side  are  known  £ priori: 

L/2  w/2 


Jx  K~+(x,y|x',y')  dx'  dy' 


0 0 


= 7 I d"  jn_1  [sin(n  + 1)4>  J , (-sp / c)  - sin(n  - l)<f>  J ,(-sp/c)]  t/a) 
y m n n-x  n i 

n=0 
n even 


5.  Baum,  C.  E. , "Interaction  of  Electromagnetic  Fields  with  any  Object 
which  has  an  Elec  romagnetic  Symmetry  Plane,"  Interaction  Note  63, 
Air  Force  Weapons  Laboratory,  Kirtland  AFB,  NM,  March  1971. 
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iMBWi -Twin iinnriT  nr  vi-r’m  - 


and 


L/2  v/ 2 


f J J ^(x.ylx'.y1)  dx'  dy' 


0 0 

00 

TTC  r .n+1 


— I jn  dx  [cos(n  + 1H  Jn+1(-up/c>  + un_^cos(n  - 1)<J  J^C-sp/c)]  (7b) 
n=0 


n even 
where 


K (x,y|x',y')  = K(x,y|x',y')  - K(x, y | -x ' ,y ' ) 

+ K(x,y|x\-y')  - K(x, y | -x ’ ,-y ’ ) 


and 


K (x , y | x ' , y ' ) = KCx.yjx'.y')  + K(x,y | -x' ,y ' ) 

- K(x,y | x ' , -y ')  - K(x,y|-x',-y') 
For  subsequent  reference 

K^Cx.ylx'  ,y')  = K(x,y|>’,y')  + K(x,y|  -x'  ,y ') 

+ K(x,y|x\-y')  + K(x,y  | -x ' ,-y ') 

and 


(8a) 


(8o) 


(8c) 


K (x , y | x ' , y ' ) = K(x,y|x',y')  - K(x,y| -x' ,y') 

- K(x,y|x' ,-y’)  + K(x,y | -x' , -y ' ) (8d) 

are  defined  ns  well.  Equations  (7)  are  enforced  for  z = 0,  x c (0,l/2) 
and  y e (0, w/2) . 

Table  1 summarizes  the  four  symmetry  cases  which  are  derived  in  the 
foregoing  discussion.  By  means  of  this  table, four  integral  equation  pairs 
can  be  constructed  in  the  spirit  of  (7)  by  replacing  the  kernels  in  (7)  with 
the  appropriate  kernels  from  the  table  and  retaining  only  the  non-vanishing 
terms  in  the  series  expansion. 

Figure  2 depicts  qualitatively  the  respective  modal  current  distributions 
for  the  lowest  frequency  natural  raasonance  exhibiting  each  symmetry. 
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SECTION  IV 


THE  NUMERICAL  MODEL 


The  integral  equation  pair  of  the  form  (7)  for  each  of  the  four  sym- 
metry cases  can  be  discretized  by  the  method  of  moments.  In  the  work  reported 
here,  two-dimensional,  subsectionally  constant  expansion  functions  were  used 
with  collocation  testing.  The  zoning  scheme  is  represented  in  Figure  3. 

The  unknown  currents  and  were  expanded  in  piecewise  constant 
functions  as  in  (ref.  3)  with  subsectioning  of  the  form  given  in  Figure  3. 
Notice  that  half-width  patches  are  used  at  the  edges  of  the  plate  so  that 
match  points  lie  precisely  on  the  edge  of  the  plate.  The  half-width  pulse 
has  proved  useful  in  realizing  the  actual  electrical  size  of  a body  in  one- 
dimensional problems  (ref.  6)*  Some  numerical  experimentation  was  also  done 
with  full-sized  pulses  on  the  edges  and  comparative  results  are  reported  in 
a later  section.  Some  difficulties  occur  in  definition  of  the  edge  of  the 
plate  in  the  present  formulation  because  of  the  presence  of  two  current 
components  which  have  the  asymptotic  behavior  given  in  (3) . This 

difficulty  is  discussed  in  a later  section. 

The  boundary  condition  J ■ 0 must  be  enforced  on  selected  patches 
at  the  edge  of  the  place  as  discussed  in  (ref.  3).  Concomitantly,  only  as 
many  d~'s  are  retained  in  the  right-hand  side  summation  in  (7)  as 
there  are  current  values  preassigned  to  zero.  The  shaded  patches  in  Figure  3 
indicate  the  selection  of  patches  where  a current  component  is  preassigned  a 
zero  value.  At  the  corner  patch,  both  components  are  preassigned  zero  values. 


6.  Butler,  C.  M. , "Integral  Equation  Solution  Methods,"  in  "Wire  Antennas 
and  Scatterers,"  Short  Course  Notes,  University  of  Mississippi, 

April  1972.  (See  also  IEEE  Trans. , v.  AP-20,  pp.  731-736,  1972.) 
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By  assigning  one  match  point  per  expansion  patch  and  by  retaining  one 
series  expansion  term  for  each  current  '.'alue  preasalgned  in  each  of  the  two 
integral  equations,  a consistent  (i.e.  square)  system  of  linear  equations 
results.  The  truncated  summation  is  taken  to  the  left-hand  side  so  that  a 
homogeneous  system  results.  The  matrix  organization  used  to  represent  these 
equations  is  given  in  Figure  h.  Table  2 defines  the  computer  variables 
noted  in  Figure  4,  primarily  for  reference  purposes  in  the  next  section. 

The  matrix  that  results  is  a function  of  the  complex  frequency  s.  A 
natural  resonance  occurs  when  s has  a value  such  that  the  matrix  has  a zero 
determinant;  hence,  the  homogeneous  system  of  equations  has  a non-trivial 
solution.  The  next  section  explores  some  algorithmic  considerations  in  the 
efficient  generation  and  manipulation  of  the  matrix. 
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Table  2 

MATRIX  FORMUIATION  PARAMETERS 


Nil 

NI2  « Nil 
NPREJ 


Number  of  match  points  on  the  zoned  quadrant  of  the  plate. 


Number  of  patches  along  the  |x|  = L/2  edge  where  J 


is  preassigned  to  zero. 


NPREI 


Number  of  patches  along  the  |y| 
is  preassigned  to  zero. 


w/2  edge  where  J 


NJ1  - NI1-NPREJ 
NJ2  = NI2-NPREI 


Number  of  unknown  current  values  in  expansion. 


NJ3 


Number  of  unknown  current  values  in  expansion. 


NPREf  31  NPREI-NPRKJ  Number  of  preassigned  current  values  (Also  the  number  of 

coefficients  retained  in  summation). 
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ALGORITHMIC  CONSIDERATIONS  IN 
EVALUATING  THE  SYSTEM  DETERMINANT 

Some  considerations  taken  into  account  in  generating  the  system  matrix 
and  evaluating  its  determinant  efficiently  are  discussed  in  this  section. 
Since  these  two  operations  must  be  repeatedly  carried  out  for  many  values  of 
s in  the  course  of  determining  the  natural  frequencies  of  the  plate,  it  is 
essential  that  clean,  efficient  computer  programming  and  coding  be  used  so 
that  execution  of  the  program  will  be  affordable.  The  volume  of  code  in  the 
algorithms  is  consistently  compromised  toward  a larger  size  in  order  to  meet 
the  following  two  time-efficient  objectives: 

1.  Avoidance  of  calculating  the  same  quantity  twice;  and 

2.  Avoidance  of  logical  decisions,  particularly  those  which  might 
be  imbedded  in  loops. 

The  program  is  discussed  in  the  context  of  the  following  major  segments: 

1.  Computation  of  an  "interaction  matrix"; 

2.  Construction  of  the  non-zero  submatrices  of  the  system  matrix 
from  the  interaction  matrix; 

3.  Computation  of  the  series  terms'  submatrix;  and 

4.  Determinant  evaluation. 

The  major  contribution  to  the  elimination  of  redundant  calculations  is 
the  one-time  computation  of  an  "interaction  matrix"  which  is  made  up  of  the 
individual  kernel  integral  terms  from  (2)  for  all  argument  combinations 
which  occur  in  the  computation.  The  subsequent  program  step  then  picks,  by 
subscript,  entries  from  this  matrix  and  constructs  the  appropriate  kernel 
from  one  of  equations  (8)  according  to  the  symmetry  conditions  being  solved. 
This  procedure  can  be  viewed  in  terms  of  the  layout  given  in  Figure  5a.  The 
terms  in  the  interaction  matrix  are  those  evaluated  for  the  match-point  as 
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shown  in  the  lower  left  with  the  source  patches  indexed  over  the  entire 
plate  to  generate  the  matrix.  Thus,  all  geometric  relationships  which  occur 
in  the  kernel  terms  are  encompassed  in  the  calculation.  Note  that  all  source 
patches  are  full  patches  for  this  calculaMon.  The  effect  of  half  patches 
at  the  edges  is  accounted  for  by  weighting  by  a factor  of  1/2  the  edge  con- 
tributions. The  kernel  integral  appropriate  to  the  symmetry  is  constructed 
by  summing  with  correct  signs  the  appropriate  elements  from  the  matrix. 

Figure  5b  gives  an  example  of  the  four  source  patches  entering  into  one 
kernel  integral. 

Differing  degrees  of  sophistication  are  required  in  the  calculation  of 
the  interaction  terms  depending  on  the  spacing  of  the  patches  for  which  an 
interaction  is  being  calculated.  For  the  self  patch,  i.e.,  the  patch  in  which 
the  match-point  resides,  the  integration  of  the  kernel  must  be  performed 
analytically  because  of  the  integrable  singularity  in  the  kernel  there. 
Appendix  A gives  a series  approximation  to  this  integral.  The  result  in 
Appendix  A is  evaluated  directly  in  the  program.  For  the  patches  adjacent  to 
the  patch  containing  the  match  point,  the  kernel  is  a rapidly  varying  but 
well-behaved  function.  The  integration  over  these  patches  is  evaluated 
numerically  by  a polynomial  approximation.  For  patches  further  separated, 
the  kernel  is  slowly  varying  and  the  integral  is  evaluated  approximately  as 
the  product  of  the  value  of  the  kernel  at  the  center  of  the  patch  and  the 
area  of  the  patch. 

Some  minor  time  economy  is  achieved  in  the  matrix  filling  algorithm, 
which  is  a four-dimensional  loop.  The  economy  is  found  in  the  form  of 
decision-free  indexing,  that  is,  the  source  contributions  from  interior 
patches,  from  |x|  = L/2  edge  patches,  from  |y|  = w/2  edge  patches,  and  from 
corners  take  on  different  forms.  Rather  than  index  over  all  source  patches 
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with  logical  decisions  implemented  to  discriminate  among  the  four  cases 
above,  four  different  loops  are  used. 

The  computation  of  the  series  term  submatrix  is  relatively  straight- 
forward. Because  the  Bessel-trigonometric  products  appear  in  two  terms  each, 
they  are  all  precalculated  and  stored  in  a vector.  The  individual  terms  are 
then  constructed  from  them. 

The  determinant  evaluation  profits  significantly  from  an  e;rploitation 
of  the  sparceness  of  the  matrix.  Either  of  two  approaches  may  be  taken  to 
the  sparce  matrix  manipulations.  One  is  to  separate  the  matrix  algebraically 
and  calculate  an  inverse  as  a composite  of  inverses  of  terms  involving  the 
submatrices.  The  alternative  approach  is  to  attack  t».a  matrix  directly  with 
Gaussian  elimination,  and  to  exploit  the  sparceness  directly  in  the  algorithm. 
The  latter  approach  was  chosen  for  the  present  purpose  because  it  is  judged 
to  be  slightly  faster  computationally  and  because  in  order  to  determine 
natural  mode  solutions  for  the  SEM  formulation,  the  homogeneous  system  of 
equations  occurring  at  a pole  must  be  backsolved.  The  algorithm  resulting 
from  the  second  approach  is  described  in  Appendix  B. 

The  determinant  evaluation  routine  itself  appears  in  Appendix  C as 


the  function  routine  CPLATE. 
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SECTION  VI 

NUMERICAL  CHECKS  ON  THE  ACCURACY  OF  THE  POLES 

The  results  of  some  numerical  checks  on  the  accuracy  of  the  pole  location 
as  determined  from  the  numerical  model  described  in  Sections  II  through  V are 
reported.  It  is  shown  that  the  model  predicts  well  the  poles  for  narrow 
strips  possessing  essentially  thin  scatterer  resonance  properties.  Difficul- 
ties occur,  however,  in  obtaining  self-consistent  results  under  different  zone 
sizes  for  plates  with  larger  aspect  ratios.  It  is  conjectured  that  the 
difficulty  occurs  because  the  subsectlonally  constant  current  representation 
is  inadequate  to  represent  the  correct  edge  behavior  for  the  currents- 
particularly  the  singular  behavior  for  the  parallel  component.  The  rationale 
behind  this  conjecture  is  discussed. 

Initial  tests  on  the  accuracy  of  the  model  were  made  for  a rectangular 
strip  with  a shape  ratio  w/L  - 1/10.  Such  a strip  has  an  approximate  equiv- 
alent dipole  whose  diameter-to-length  ratio  is  1/IOit. 

Figure  6 gives  the  results  of  pole  determinations  for  the  first  two 
poles  for  various  numbers  of  pulses  in  the  expansion  of  the  current  and  for 
two  different  treatments  of  the  edge  pulse.  The  strip  was  zoned  with  one 
pulse  across  a quadrant.  The  numbers  indicated  in  the  figure  are  the  numbers 
of  pulses  along  the  longitudinal  direction  of  a quadrant.  The  "half-pulse" 
results  are  those  obtained  by  the  zone  scheme  described  in  Section  IV  where 
half-width  pulses  are  placed  at  the  edge  so  that  match  points  fall  at  the 
edge.  The  "full-pulse"  results  are  those  obtained  by  zoning  the  plate  with 
equal-sized  pulses  over  the  entire  quadrant.  In  the  latter  case  air 
a,  posteriori  adjustment  is  made  in  the  data  correcting  the  length  of  the 
strip  such  that  the  end  of  the  corrected  strip  lies  at  the  end  match-point. 
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x Half  patches  at  edge 
+ Full  patches  at  edge 
with  size  correction 

• Equivalent  Cylinder 


Figure  6.  Calculated  Pole  Locations  for  Thin-Strip  for  Varying  Numbers 
of  Zones  in  the  x-Direction  and  Different  Edge  Treatments 
(Cylinder  Results  from  Ref.  6) 


It  is  seen  that  the  differences  are  small  both  for  varying  order  and 
Increasing  pulse  density.  The  NX  ■ 6 results  for  the  second  pole  show  some 
departure  from  the  trend  established  by  the  results  for  NX  » A and  NX  ■ 5. 

This  is  attributable  to  the  fact  that  the  matrix  is  on  the  brink  of  numerical 
instability  for  NX  ■ 6.  The  results  for  NX  ■ 7,  which  are  not  shown,  are 
observed  to  be  meaningless  because  of  the  instability  manifested. 

For  comparison  purposes,  the  first  two  poles  for  an  equivalent  cylinder 
(one  whose  circumference  equals  the  strip  width)  are  given  as  found  in 
ref.  7.  These  results  are  judged  reliable  inasmuch  as  they  have  been  cross- 

checked among  several  integral  equation  formulations  and  for  several  method- 
of-moments  schemes.  The  equivalent  radius  taken  is,  of  course,  an  approx- 
imation. It  is  seen  that  the  half-pulse  solutions  compare  slightly  more 
favorably  with  the  cylinder  results.  Because  of  this,  and  moreover,  because 
the  a posteriori  data  adjustment  is  avoided  with  the  half-pulse  scheme,  it 
was  used  consistently  in  the  remaining  solutions. 

The  stable  convergence  properties  of  the  numerical  model  exhibited  for 
the  thin-strip  solution  are  not  manifested  for  higher  aspect  ratios.  The 
reason  for  the  difference  is  that  the  strip  is  essentially  a one-dimensional 
problem  and  the  transverse  (y-directed)  component  of  current  has  little 
effect  on  the  dominant  longitudinal  current.  For  wider  structures  the  coupling 
is  significant  and  inadequacies  in  the  modeling  of  the  singularities  of  the 
currents  produce  inaccuracies  which  are  highly  sensitive  to  zoning. 

Figure  7 shows  the  results  obtained  for  a pole  trajectory  as  a function 
of  the  shape  factor  w/L  where  the  zoning  in  the  y-direction  was  adjusted 

7.  Umashankar,  K.  R.,  "The  Calculation  of  Electromagnetic  Transient  Currents 
on  Thin  Perfectly  Conducting  Bodies  Using  the  Singularily  Expansion 
Method,"  Ph.D.  Thesis,  University  of  Mississippi,  pp.  33-34,  August  1974, 
(See  also  Tesche,  F.  M.,  IEEE  Trans. , Vol.  AP-21,  No.  1,  pp.  53-62>  1972.) 
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according  to  the  value  of  w.  It  is  evident  that  the  solutions  are  unstable 
with  respect  to  the  zoning  on  the  plate.  Attempts  to  increase  the  number  of 
zones  significantly  to  improve  upon  the  situation  resulted  in  numerical 
Instabilities  in  the  matrix  which  cause  the  pole  search  iteration  to  fail. 

The  reason  for  the  difficulty  manifested  in  Figure  6 is  believed  to  lie 
in  the  way  that  the  edge  of  the  plate  is  defined  with  the  piecewise  constant 
current  expansion.  Consider  the  characteristics  of  the  two  current  components 
along  a line  traversing  the  plate  in  the  y-directicn  as  shown  in  Figure  8. 

The  correct  edge  behavior  at  |y|  ■*  w/2  is  that  given  in  equations  (3).  The 
zoning  scheme,  however,  forces  Jx(x,±w/2)  to  take  a finite  value.  The  current 
extrapolates  to  a singular  point  for  some  y > w/2,  i.e.,  the  numerical  model 
represents  a plate  whose  width  is  greater  than  w. 

If  the  number  of  transverse  zones  is  increased  as  indicated  by  the 
dashed  curve  in  Figure  8,  the  steepness  of  the  edge  behavior  of  is 
increased,  and  the  extrapolation  is  characteristic  of  a narrower  plate  as 
compared  to  the  first  case.  This  narrowing  of  the  effective  width  in  the 
model  is  reflected  in  an  increased  Q (resonance  quality  factor)  as  the  jumps 
in  Figure  7 indicate. 

One  is  tempted  to  conclude  that  a full-width  pulse  at  the  edge  is  a 
cure  for  the  problem  since  the  point  at  which  the  pulse  current  is  defined 
is  shifted  relative  to  the  edge  as  zoning  is  changed  with  full-width  pulses. 
The  effect  of  this  procedure  is  to  transfer  the  problem  from  component  of 
current  whqse  behavior  is  singular  at  the  edge  to  the  component  which  goes 
to  zero.  With  full  pulses  at  the  edges,  the  normal  component  of  current 
would  go  to  zero  interior  to  the  plate  rather  than  at  the  edge  as  it  properly 
should. 

An  approach  which  is  potentially  a remedy  for  this  difficulty  is 
discussed  in  the  conclusions. 
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SECTION  VII 

POLE  TRAJECTORIES  AS  A FUNCTION  OF  SHAPE  RATIO 

Figure  9 gives  the  results  obtained  for  pole  location  for  the  lowest 
order  pole  of  each  of  the  symmetries  as  a function  of  w/L.  Clearly,  as  the 
previous  section  indicates,  the  results  cannot  be  taken  as  the  correct 
results  for  the  plate.  However,  the  zoning  was  fixed  for  all  w/L  and  the 
results  are  expected  to  reflect  the  proper  trends  for  these  trajectories. 

The  ++  and  +-  modes  are  in  essence  dipole  modes,  and  their  poles  show 
the  general  lowering  of  0 as  w/L  increases.  (The  ++  indicates  that  the 
component  is  symmetric  both  with  respect  to  the  x and  y axes  - see  Table  I.) 
The  — .node  can  be  thought  of  as  a dipole  mode  in  the  transverse  direction. 

As  a result  if  shows  a strong  frequency  dependence  on  the  transverse  dimension 
w.  When  w/L  = 1,  the  — mode  is  identical  to  the  ++  mode  rotated  spatially 
90  degrees.  Consequently,  the  two  trajectories  coalesce  as  w/L  1,  when 
the  zoning  is  5x5  so  as  to  preserve  symmetry  in  the  numerical  mode.  The 
failure  of  the  5x3  zone  case  is  due  to  the  reasons  outlined  in  the  previous 
section. 
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SECTION  VIII 


CONCLUSIONS 

The  application  of  SEM  to  the  equivalent  problems  of  the  perfectly 
conducting  rectangular  plate  and  the  rectangular  aperture  in  a perfectly 
conducting  screen  has  been  conducted  with  partial  success.  The  applicability 
of  SEM  and  the  computational  feasibility  of  determining  SEM  quantities  are 
demonstrated.  It  is  further  demonstrated  that  by  careful  program  construc- 
tion, the  computational  costs  of  numerical  treatment  of  two-dimensional 
problems  can  be  made  quite  reasonable.  The  cost  of  generating  a matrix  and 
calculating  its  determinant  by  the  methods  discussed  herein  is  J.ess  than 
10  cents  for  each  value  of  8. 

Difficulties  arise  in  the  subsect ionally  constant  current  zoning  be- 
cause of  a failure  to  properly  model  the  edge  conditions.  Whereas  Rahmat-Samii 
and  Mittra  (ref.  3)  obtained  good  radar  cross-section  results  with  such  a 
zoning  scheme,  the  pole  locations  are  highly  sensitive  to  the  zoning.  Radar 
cross-section  is  a far-field  quantity,  and  the  integrated  effects  of  the 
errors  are  small  there.  The  pole  locations,  on  the  other  hand,  depend  strongly 
on  the  dimensions  of  the  structure,  and  it  is  evident  that  the  plate  size 
must  be  brought  to  bear  in  a precise  fashion  for  the  pole  locations  to  be 
correct. 

The  edge  problem  can  be  remedied  by  imposing  the  edge  conditions 
(3)  directly  in  the  basis  set  elements  for  edge  zones. 

Davis  has  done  this  successfully  for  the  circumferential  component  of  current 
on  a thick  cylindrical  scatterer  (ref.  8).  The  circumferential  current 

8.  Davis,  W.  A.,  "Numerical  Solutions  to  the  Problem  of  Electromagnetic 

Radiation  and  Scattering  by  a Finite  Cylinder,"  Ph.D.  Thesis,  University 
of  Illinois,  1974. 


component  is  singular  at  the  ends  of  the  cylinder.  The  introduction  of  the 
singular  basis  element  will  produce  a significant  complication  to  the 
problem  in  that  a second  singularity,  that  of  the  current,  must  be  inte- 
grated analytically  in  order  to  implement  the  model  with  edge  conditions 
imposed.  An  independent  check  on  program  accuracy  is  dictated  for  a 
problem  of  this  scope  before  proceeding  with  the  edge  condition  approach. 


APPENDIX  A 


THE  SELF-PATCH  INTEGRATION 

The  term  for  the  interaction  matrix  for  IDIF  ■ JDIF  * 0,  i.e.,  where  the 
match  point  lies  at  the  center  of  the  source  patch,  can  be  written 

Ax/2  Ay/2 

I K(0,0 | x' ,y ' ) dx’  dy'  (Al) 

0 0 

This  presumes  a unit  amplitude  expansion  pulse  over  the  patch  whose  dimensions 
are  Ax  and  Ay.  The  symmetry  of  the  kernel  with  respect  to  both  x'  and  y'  is 
employed  in  the  forming  of  (Al).  This  integral  can  be  transformed  to  polar 
coordinates  as 


1-4 
s H 


. -1  Ay  Ax 

tan  , 7 

j Ax  2c^s<ji 

<t>=0  p-0 

Ax 


exp[-sp/c]  dp  d<j> 


2cop<J> 

<t>=tan  p=0 

Ax 


exp[-sp/c)  dp  d4> 


-1  A£ 


itan  t* *• 
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[exp(-sAx  sec  4>/2c)  - 1]  d<J> 


$=0 


tt/2 


+ I [exp(-sAy  esc  4>/2c)  - 1]  d<j> 

V1  & 


<J>=tan 


Ax 


(A2) 


If  the  exponential  functions  in  the  integrand  are  then  expanded  in  McLauren 
series,  the  resulting  terms  of  powers  of  secants  and  cosecants  possess  tabu- 
lated integrals.  The  result  for  three  terms  retained  in  the  series  is 
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APPENDIX  B 


THE  SPARSE  MATRIX  ALGORITHMS 


1.  Introduction 


This  Appendix  describes  the  algorithmic  approach  to  minimize  the 


computation  time  involved  in  Gaussian  elimination  triangularization  of 


systems  of  matrix  equations  which  are  "sparsely  coupled."  The  term 


"sparsely  coupled"  as  applied  in  this  Appendix  implies  the  matrix  equation 


form  given  in  (Bl). 


[M]  [X] 


It  is  clear  that  in  this  form  the  sole  coupling  between  the  "upper"  and 


"lower"  systems  of  equations  is  contained  in  the  matrix  M,,.  Generally,  the 


number  of  columns  in  M2  is  small  compared  with  the  order  of  the  overall 


system. 


An  algebraic  approach  to  exploiting  the  sparceness  of  (Bl)  results  in 


solutions  which  are  given  in  terms  of  the  several  submatrices  and  their 


inverses.  (See,  for  example,  ref.  9.)  It  is  well-known,  however,  that  it  is 


sufficient  for  the  purposes  of  determinant  calculation  and  equation  solution 


to  triangularize  the  composite  matrix  in  (Bl) . The  triangularization  process 


involves  fewer  operations  than  the  diagonalization  necessary  for  the  calcu- 


lation of  an  inverse. 


9.  Dunaway,  0.  C.,  "A  Numerical  Solution  for  the  Distribution  of  Time- 
Harmonic  Electromagnetic  Fields  in  an  Arbitrary  Shaped  Aperture  in  a 
Ground  Screen,"  M.S.  Thesis,  University  of  Mississippi,  1974. 
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|0*I1  x NJ1) J 


(Nil  x NJ2) 


! o 

CMAT2 

(NI2  x NJ1) 

(NI2  x NJ2)‘ 

a.  _ 

CMAT3 

(NI3  x NJ3) 


NI3  - Nil  + NI2 
NJ3  = NI3  - NJ1  - NJ2 
(a) 


NI4  = MAX  (Nil  - NJ1,0) 
(b) 


Figure  Bl.  Submatrix  Organization  for  the  Sparse  Matrix  Algorithms,  a) 
the  Original  Matrix,  and  b)  Triangularized  Form  with  the 
Generated  CMAT4 
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This  Appendix  describes  an  algorithmic  exploitation  of  the  sparseness  of 
the  composite  matrix  in  (Bl)  • That  is,  a numerical  process  is  given  whereby 
only  the  non-zero  subelements  are  stored  and  operated  on,  with  the  compu- 
tational operations  being  those  which  render  the  composite  matrix  M upper 
triangular.  The  determinant  of  the  composite  matrix  results  directly  from 
this  triangularization.  A solution  for  X in  (Bl)  requires  a backsolving 
process  involving  the  triangularized  form  of  M and  a vector  resulting  from 
applying  the  elimination  operations  to  the  vector  B . Routines  to  perform 
these  operations  are  given  as  well. 

Appendix  C gives  listings  of  the  routines  built  on  this  algorithm.  The 
triangularization  routine  is  named  SPRHOM.  The  backsolving  procedure  is 
performed  by  the  entry  HOMSLV  to  the  routine  SPRSLV.  (The  name  entry  SPRSLV 
backsolves  an  inhomogeneous  system  and  is  not  us  -A  for  present  purposes.) 

2.  The  Algorithm 

The  routine  SPRHOM  is  simply  an  implementation  of  the  Gaussian  elim- 
ination procedure  with  maximum  pivot  selection  applied  to  the  composite 
matrix  M in  (Bl)  viewed  as  a whole.  The  sparseness  of  M is  exploited  by 
storing  only  the  non-zero  submatrices  in  (Bl)  and  skipping  the  operations 
involving  zero  elements.  The  result  is  a substantial  saving  in  both  storage 
and  computation  time. 

The  forms  of  the  input  and  of  the  end  product  for  the  triangularization 
are  given  in  Figure  (Bl)  with  the  sizes  of  the  respective  submatrices  defined. 
It  is  recalled  that  the  fundamental  process  in  the  Gaussian  elimination 
procedure  is  the  elimination  of  all  or  part  of  the  elements  of  a column  of 
a matrix  with  respect  to  a pivot  element,  commonly  the  element  of  greatest 
magnitude  in  the  column.  That  is,  for  a column  under  process,  the  row 
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containing  the  main  diagonal  element  of  the  matrix  which  falls  in  that 
column.  All  or  part  of  the  elements  not  on  the  main  diagonal  are  "eliminated" 
or  made  zero  by  subtraction  of  some  multiple  of  the  row  containing  the  col- 
umn maximum.  In  the  triangularization  procedure,  the  part  of  the  column  com- 
prising elements  lying  below  the  main  diagonal  after  row  exchange  are  elim- 
inated. If  the  matrix  is  a part  of  a system  of  equations  with  non-zero  right- 
hand  side,  the  row  operations  of  exchange  and  subtraction  of  a constant 
multiple  of  another  row  must  be  performed  on  the  corresponding  elements  of 
the  right-hand  side  vector  as  well. 

The  algorithm  of  the  routine  SPRHOM  operates  according  to  the  Gaussian 
elimination  procedure  described  above.  However,  the  three  submatrices  CMAT1, 
CMAT2,  and  CMAT3  are  stored  individually.  In  addition,  the  routine  generates 
a submatrix  CMAT4  in  the  course  of  selecting  pivots  for  the  columns  contained 
in  CMAT2.  Further,  the  "elimination"  of  elements  of  submatrices  that  are 
zero  a priori  is  skipped.  The  result  Is  significant  storage  and  time  economy. 

In  order  to  minimize  logic  decisions  in  the  routine,  it  is  organized  to 
operate  sequentially  through  the  partitioned  matrix.  The  steps  are  as  follows 
(it  is  convenient  to  follow  the  thinking  of  these  steps  by  tracing  the  loca- 
tion diagonal  of  the  composite  with  the  aid  of  Table  Bl)  : 

a.  Perform  conventional  Gaussian  elimination  to  zero  the  elements 
CMAT1(I,J)  for  I > J,  i.e.,  the  elements  below  the  main  diag- 
onal of  M.  Choose  maximum  pivot  elements  in  conventional 
fashion.  Carry  row  operations  into  CMAT3. 

b.  Create  CMAT4  by  row  swapping  with  CMAT2  so  as  to  choose 
maximum  pivot  elements.  Perform  elimination  to  zero  CMAT4 
elements  for  I > J and  the  entire  column  of  CMAT2.  The 
number  of  rows  created  in  CMAT2  is  NI4  = Nil  - NJl,  the 
amount  by  which  CMAT1  is  over-square.  Carry  row  opera- 
tions across  into  CMAT3. 
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c.  Choose  maximum  pivot  rows  in  columns  of  CMAT2  with 
J > NI4  and  swap  with  rows  given  by  I ■ J - NI4 

(the  rows  containing  the  Jth  column  diagonal  element  of 
the  composite).  Conduct  elimination  to  zero  elements 
with  I > J + NI4.  Carry  row  operations  into  CMAT3. 

d.  Conduct  conventional  pivot  selection  and  elimination 
on  CMAT3  to  zero  elements  CMAT3(I,  J)  with 

I > J + NJ1  + NJ2. 

In  each  column  elimination  operation,  the  pivot  value  is  multiplied  into  a 
product  accumulator  to  produce  a value  for  the  determinant  of  the  composite 
matrix.  The  sign  of  this  product  is  changed  at  each  row  swap  in  accord  with 
the  rates  of  matrix  algebra  row  operations. 

The  backsolving  routine  SPRSLV  with  its  entry  HOMSLV  operate  in  a 
straightforward  manner  on  the  submatrices  as  reduced  by  SPRHOM. 

Details  are  omitted  here,  but  the  routines  may  be  easily  followed  in  a row- 
by-row  flow  from  the  bottom  of  the  composite  matrix,  if  one  keeps  in  mind 
the  row  index  relations  of  column  4 of  Table  Bl.  The  entry  HOMSLV  assumes  a 
zero  determinant  value  resulted  (approximately)  from  SPRHOM  and  the  last 
element  of  the  solution  vector  is  picked  as  unity.  (The  zero  determinant 
results  from  a zero  falling  at  the  last  diagonal  location  in  maximum  pivoting 
Gaussian  elimination.)  The  remainder  of  the  homogeneous  solution  vector  is 
backsolved  conventionally  relative  to  this  last  element.  The  vector  is  then 
renormalized  so  that  its  maximum  element  is  unity. 


Table  Bl 
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PRIMARY  INDEXING  QUANTITIES  IN  THE  ALGORITHM 


Submatrix 


Size  of 
Submatrix 


Indices  of  Main 
Diag.  of  Compos. 


Relative  Row  Index 
of  CMAT3  and  CRHS 


CMAT1 

Nil  x NI2 

(J,J) 

CMAT4 

Nil  - NJ1  x NJ2 
(can  be  null) 

(J,J) 

CMAT2 

NI2  x NJ2 

(J  - (Nil  - NJ1),  J) 

CMAT3 

Nil  + NI2  x 
Nil  + NI2  - 
NJ1  - NJ2 

(J  + NJ1  + NJ2,  J) 

13  = I + NJ1 


1.  Quantities  given  in  terms  of  input  parms.  to  the  routine.  Related 
internal  quantities  are  given  in  Figure  Bl. 

2.  Relative  to  I,  the  row  index  of  the  submatrix  in  question. 
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APPENDIX  C 


PROGRAM  LISTINGS 


All  code  compileable  on  IBM  OS/360  and  OS/370  FORTRAN  levels  G or  H. 


The  routine  ZANLYT  and  its  service  routine  UERTST  is  taken  from  the  program 


library  FORTUOI  made  available  by  the  Computer  Services  Office,  University  of 


Illinois  at  Ur b ana-Champaign,  Urbana,  Illinois  61801.  The  routines  BSUZ 


and  BSCJZ  are  taken  from  the  International  Mathematical  and  Statistical 


Library  (IMSL) . They  may  not  be  reproduced  apart  from  this  application 


program  package.  The  IMSL  library  is  available  by  subscription  from  IMSL, 


Inc.,  6100  Hillcroft,  Suite  510,  Houston,  Texas  77036. 
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■■  Iff 

1 


******** 


• PHLF  RFARCH  oprGOAM  Fp»  S F M FTP m'JlATTPM  OF  THIK-PLATF  SCATTFBFR 
BY  L 4 PFARSPN  0/74 

I MPLI CT T RFAL*8( A ,P,P-H,r-ZI ,PPMOLFX*16 ( F ) 

C?mmOS  /GFOM / XSYMtYSYw»WfNXf0Y,I0PFAS( JO  ) »JPPFAF(10)fNpoF!»MPOFJ 
TNT  FGFP  MFS ( 4 • 2 ) / 1 SYMWt , t cTb  T • , *r  i , i • , «A NTI  ' S VMM  • t ' E*  R I * * ' r ' / 
DATA  C /(  3. 008, 0. OP)  /,0|.  US/ •«■•/,  PI/3.141592653509793/ 

DATA  MX/’X'/t HY / • Y • / 

EXTFOMAL  r“PLATF 
DIMENSION  :x(20»f INFFP<20) 
logical  LAJTO 

100  RF  AD(  5*1  ,FN0*999)  X SYM  , Y SYM,  NX  , t'Y  , WO,  W£  ,W«  , r$(JNrp  , I.  AUTC 
1 FOR  MAT  ( ‘»Al,2X,2n,5Fin.4,TBft,L’  ) 

I MX  * 1 
IMY*l 

»F(XSYM.NF.PL'JS  ) TMX  = ? 

TF( YFYM.NC.PLUP ) T UY  =2 
\JW*(WM-WP)/WS 
! F ( AM  ,0T  .0  ) OP  T°  10^ 

NW=-NW 
WS=-WS 

105  !F(WS*NW.LT.WM-Wr»  ) NW«NW*1 
DP  20*  I W=1 
W*WP*( TW-1 )*W$ 

IFf  .NPT.LAjrn)  r,n  -n  1<vp 

S*IP  P/ST  AUTO  ZnN INO 

P*UT I N F TO  DFTCPMTNF  NH  op  FXPANSI*N  Pl'lSES  B^SFO  PN  CI.FCTBIr  AL 
SIZE  hf  PL'tF 

TCS ‘»,WV*«  1.8850 10/D  ARSIDt* *01  * SUNPP  ) ) 

////////////// 

NPPWV  La6 

////////////// 

FLFNX*1/TFSTWV 

NX  = I 01  NT ( FL  FN  X*NP°  WVL  ) 

tF{0F1.0AT(NX)  .IT.FLFNX*NP°WVI  ) MX=NX  + 1 
r.LFMY=W/TESTWV 
'vY=IOTNT(FLENY*MPPWVL  I 
V F(DFLr'AT(NY)  .IT.  FLF^y*NPPWVU  NY  = MY  + 1 
NX  = mIN0(NX,5) 

NY=mtnO(NY,51 

r 

* RFP,IN  SFTUP  F no  ALTFPNATf  Epr.F  P/TCH  °o F A S S I ONM=N~ 

C 

140  NPPET  =( NX+-7  ) /3 
NPR  FJ  = (N Y+2  • / 3 

T F ( nx-2*NPP  F I *2  .L  c . 1 . am 0 .MOO  c i . GT  . 1 ) N PF F T =0PB  B I - 1 
Ic(MY-2*NPRej«-2.LF.l.»ND.MORPJ.GT.l  ) NPRFJ=NPP  FJ-I 
OP  110  T«1.NPPFI 
TPPFAMNPRFIM-U  =NX-3*T+3 
HO  "ONTINUF 

OP  120  J-l.NPRFJ 
JRRFASINPRF J+l-J) = N Y-  3*J*3 
120  CONTINUE 

C LPPA'r  T^N’S  WHE°  F X-DIPFCTFO  ruDPcN 

f T IF  £FT  TP  ZFor'  PIVFA'  ay  cyRFn] 


00010 
00020 
*0*30 
00040 
00050 
*0060 
00070 
00080 
00090 
OOIOO 
oo]  1 o 
00120 
001  3'' 
o*l  4* 
001  50 
00160 

001  70 
00180 
00190 
oo?  on 
'">0?  1 0 
00220 
O023P 
00  2 40 
00250 
00260 

002  70 
00280 
00290 
00300 
00310 
00320 
00330 
00340 

003  50 
00360 
00370 
*03  BP 
00^90 
00400 
0*410 
00420 
00430 
00440 
00450 
00460 

004  70 
00480 
0 04  9* 
00500 
00510 
*0520 
00530 
00540 
00  5 50 
00560 
nns  7 n 

0058* 
00  5 90 
00600 
006  10 
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PTS  INX.JPREAS)  AND  Y-T IB  FCTFO  BY  00620 
( I PF  EAS  * MY  l 0063P 

00640 

WPITE(6*2J  W,r<UNPP  00650 

C3RMAT<  • 1FNTEP  !tPPATTPN*  ,/t'nSHAPP  PAMP  *',F5.3,6X,  00660 

PSTAPTINO  FIFO  «',2ni2.4)  00670 

MPTTP(6,3I  00680 

p-i»M  AT  ( • O'  * 10 X*  ,rUP  SYMMFTBy*  ,6X,  • PUl  SFS' ,3X,*  PPEASS' LPP"NS')  00600 
W°  I T E ( 6 , 4 ) HX,<MESI  T,mi,T*l,4),NX,l,  tPPEfiS<J>,J*l,NPPEII  00  7 00 

pop  MAT  ( • • , A 1 ,'-0!B',5X, 4*4, !6,5X, 1013)  00710 

WPTTP(6,A)  HY, <MPS( ! ,IMY) ,T*1  ,6) ,NY,(  JPOFAS( J) , J*!  ,NPPFJ)  00  7 20 

WRITP ( 6»  5 ) 00730 

FnPMAT(»o*,llX, 'COMPLEX  FPFO' ,1 7X , • 0FT FPM I NA  NT  • ) 00760 

CX(l  IsOFlJN^P  00750 

CALL  ZAMLYTICPL ATE, 1 .0-50,4,0, 1,  l,rx, l^O, INFER , IFP  ) 00760 

W3  T Tc ( 6,6)  :x(1  J 00770 

FPPMAT{ 'OPETJFN  FP.rw  ITE°ATTnN' »/*  'OPCLE  LCC ' • N ' , 2E l 2. 4)  00780 

0 ALL  “OOF  00700 

csuN3P*rxm  oosoo 

CONTINUE  0081 P 

S3  TO  100  00820 

STtp  00833 

END  0*840 
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SUBonjTTNE  MOOt'  0085C 

T MOLT  TIT  BEAL*8<A,*»0-H,r'-Z),rnMPt  ?X*16(C>  00860 

COMMON  /MAT / CMATX<  25,?5)  ,CM»TY<25v..iu»  ,CH0«(50,10)  ,CM/T4(10,25),  0087C 

inptchs  ,noimi  , noimc i,NnjMCj,N"ep  ooaac 

COMM^l  /GEHM/  XSYM,YSYM,W,»JXtNY,JPRr#SnO),JPPEASC10),NPPEI,NPRFJ  00890 
9tMENST0N  rPRX(5,5),CPPY<5,5|  00900 

DIMENSION  C J (50 1 00910 

Ni>PE.NPREI*NPPFJ  00920 

NPO  EIM*NPRFT- 1 00930 

NPREJM-NPREJ-1  00990 

CALI  HOMSLV(CMATX,NPTCH«;,NPTrHS-NPREJ»M''IMlf  NDIMl  , 00950 

1 f MATY  * NPTC  H5  » MPTf,H  G-NPR  El  »NDI M 1 , ND I Ml , 00960 

2 '•HnMtNni‘i<‘T.NPTMrj,rMAT4tNniMCJ,NOIMltCJ»NORO)  00970 

NXM1.NX-1  00980 

NYM1«NY-1  00990 

NSUPS-0  01000 

00  470  JS=1,NY  01010 

Of  450  I S*1  * NXM 1 01020 

J=( JS-l)*NX*T $ 01030 

C°RXUS,  JS)«CJ(  J-NSlJRS)  01040 

JM*J-NSU8S  01050 

COMTINUF  01060 

J* JS*NX  01070 

TF  ( JS.NE.JPPE  AS  (N$U«SM)  ) GP  T0  460  01080 

N«UP$*M!N0(NSUMS+1,NPR1 JM|  01090 

CPRX(NX,JSI  = (0«,,0.)  01100 

GO  TO  470  OHIO 

CP»X(NX,JS) *C J(  J-NSU3S)  01120 

CONTINUE  01130 

00  500  T$«l,NX  01140 

03  500  JS*1,NYM1  01150 

J«(JS-1I*NX+1S  01160 

rPPY(T5, JS1*C  JINPTrHS-NPPFJ*J)  01170 

CONTTNUE  01180 

NP  |PS*0  01190 

On  53n  IS-lfNX  01200 

J=NYM1*NX*1 S 01210 

1 FI  TS  .NF.IPPE  AS  (NSU«S«-1 ) ) r,p  Tn  510  01220 

CPRY(TS,NY)=(0.,0.)  01230 

NSUPS«rtIN'0(  NSU8S*-1»NPPETM)  01240 

GO  T9  530  01250 

r.PR  y(  i s»  uy)  *r  j(nptchs-nppej*j-nsurj;  01260 

CONTINUE  01270 

MR !TP ( 6 , 16 ) 01280 

FnjHATJ •0*****NATURAL  MrDc***v* » ,/ t 'OX-OI  R FCTEO  CURRENT')  01290 

00  540  1 = 1,  NX  01300 

WPITF(6,17)  (CPRX(!*J) » J = 1 *NY ) 01310 

FTP  MAT ( * O'  * 51  2012.4*2X))  01320 

CPNTTNUF  01 33C 

W®T  TE (6  * 181  01340 

F3RMAT(  'OY-DIOEPTFn  'UR5FMT'J  01350 

0"  550  1=1, NX  0136^ 

WPT  TF  C 6,1  7)  (rPPY(I,J)  ,J«1  ,NY)  01370 

CONTTNUF  01380 

M°ITF(6,19)  0139" 

FORMAT!  'OHOMOGENFCUS  EXPANSION  CPEF"S'I  01400 

MPTTF!6,17)  (CJ(2*NPTCMS-WPR=  + n,I»l,NPPE)  01410 

RFTU°N  01420 

FND  01430 
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COMPLEX  FUNCTION  C»l * T *16 ( CsuNQe ) 

0=TFRHtNANT  E VALUATION  PCUTTME  FPR  HALL  FN-TY  PF  All  GMFNT  FD  MOMFN’ 
'iATRtX  Fpc  THe  THIN  Pt  AtF  SC  ATT  FP  EP 
FOP  SEN  A°PLKAtTPM« 
py  l W °FA®  SON  8/74 


IMPLir  IT  rpMPLFX*16(PJ  ,<>  F»l*8(  A ,B  , D -H,  P-Z  ) 

GOYMON  /S  CDV  XSYM,  YSYM,  W, NX, NY, I PP^AS  UD  ) , JPPFAS  ( 1 0 ) , NPRE! , MP«FJ 
C.PWM7M  /MAT/  r MATX(  25,25),CMAtY<  25,  2 5)  ,PHPM(  50,10)  *CYAT4(  10,25)  , 
lNPTHC.UDJMl  , NDI  Mp  T ,wr>T  mt  J , NPPP 


l-NV-M  jINIIl"1!.!  , I'm;  i -i  J , T 

1P4L+8  OPA^O.OTMAFG.OPPFFl’O)  ,0! MR F $ ( 20 ) , OUM]  (20)  ,0MM2  ( 20  ) , 01IM3  ( 20 
1)  ,OU”4(20) 

51  MrNSI  PN  CT'TP  (10,10)  ,r  TNTX<25)  ,C  INTY(?5  ) ,CCPSTW(  10)  ,GS  1^TM  C 10) 
j\|TrOFi?  ME5(4,2  ) / ' SYMM* , ,eT5  T ',  'f  ' • , 'A  NTI  • , • SY^M*  , • ET»  1 ' , • C • / 

5 <T*  0/(3.008,0.00)/, PL  US/' ♦,/,PT/3.141  59  2653589793/ 

kir»  " it  I _ 'i  c 


NO" Ml»?5 
NOJT 1*50 
NOTMO  J*n 
90TM=50 


fpomijl?tION  f ftijp  pput!mF? 


TMX*l 
I YY  = 1 

I F(  XSYM.MF.PLUS)  T^X=-> 

IF(Y5YH.MF.£»Lii«:  ) T»Y=? 


JV  (X/Y  )~2  TNDTATS  f NTI  SYWTTP  IC 
01STR  OF  X-DI P FPr  FO  CURF  ENT  WPT  X 
/Y  AXIS 


N°TC  H P--sMX*MY 


TFT  NO  r F CUPP  FNT  P ATCHES 


NXM1*NX-1 
MYM1  =9  Y- 1 
EPGFA  0 =0..  5 


pepprcTIONAL  W T DTH  OF  EDGE  PULSE* 


EOG?*FDGFAr*EOGFAr 


CORNER  FACT  PF 


DX=1  . / (FLPAT(  2*MX-2  )*2*Fnr.FAr  ) 

"Y=W/(  FLPAT(  2 *MY-2  ) ♦?*Fn<*,F  A*  ) 

\IXT2=MX*2 

NYT?=MY*2 

l*S  = CSUNno  /?/" 


NTTL  TZFO  LAPLAOF  VAPTAPLE 


T 9TPT5  =1 3 
nXTA'T  = 0X/12 
0YINT=r)Y/l2 


NUTR  I A'T  F G PA  F MS 


9*YMX=-(-l )**TMX 
MSymy=-(-1)**!my 


SIGNED  SYMMFTRY  TNDTT ATpO* 


NPmtt*95yMY 

MSmti 1=NSYNX*NFY9Y 

NSMIV®N$YMX 


SIGNS  CF  KERNEL  F"o  FA  OUAD'S  CON 
TP  T PUT  TPM 


NTNDX=  2 

TF(fSMIII.GT.O)  N T M OX  = 1 


NINOX  = I INDTGATFS  FVEN-FVFN  "F 
PDO-PDO  SYMMfTPY  FPP  X-DIR  CUPP 


NSCCC=1 

I F ( X SYY  .c 0.  Pt  US)  nE.fpt  = t 


01440 

014  50 
01460 
01470 
01480 
01490 

015  00 
01510 
01520 
01530 
01540 
01 5 50 
01560 
01570 
01580 
01590 
01600 
01610 
01620 
01630 
01640 
01650 
01660 
01670 
01680 
01690 
01700 
01 7 1C 
01720 
01730 
01740 
01750 

017  60 
01770 
01780 
01790 
01800 

018  IP 
018  20 
01835 
01840 
01350 
01R6P 
01870 
01880 
"1  390 
01900 
01910 

019  20 
01930 
0194P 
01950 
01960 
P197P 
01980 
01990 
02000 
02010 
02020 
0 2"  30 
02040 


■>  o n •)  ft  .1  -»  -»  «r»  v~» 


iy>pHRFniiil.i,M  i|.  . ll.HIWf 


PP 


C 

r 


r 

c 


\iscrs  » 2 i n'oic ^ te s even  symm  hrt 

Y FT R X DIR  CHRP  <1  F COSINE  EXP* 
NS  I PN  rp  HnNCENFOMS  Sni'*l) 


NPRF»NPREI»NPREJ 


TPT  Kir  CF  PPE&SSTONEO  CUPP  VHHS 


NPRFJM»NPRFJ-1 

nprfim-nprfi-i 

\(PRFPl«NPPE4-l 

FNO  OF  TNPUT/SPECIFIC«T IPN  ROUTINES 


BPUTINF  to  fill  TNTFF^CTIHN  matrix  FPPM  WHICH  MOMFnt  matrix  is 
r.ONSTPUCTFO 


OTAO,*p$nRT(nx*DX«-PY*OY) 

ALNXTM»2*OLOG ( ( 01 AG+OY )/OX  » 

9LNYTH«2*DL?G ( ( DI  AO+0X ) /QY ) 

OY8PX*0Y/0X 

0X*nY=0X/0Y 

c$ox*:s+dx 

C:0Y*CS*DY 

:<JOX2«CSDX*CSOX 

csox?*csox*csox? 

CSOX4«C$dX2*CSDX2 

CS9Y2*CSOY*CSOY 

CSDY3*CSDY*CS0Y2 

CS0Y<t*CS0Y2*CSPY2 

CrMPnNFKIT  TfCMS  Fro  sHF-PAtch  SE 

RIFS 

CXTFBM»-0.'>nP$CSOX*fiLNXTM+''.5nO*CSDX2*DYPDX-C'DX3*(  0XPPY*D  JAG/ll  2* 
10Y1  + ALNXTM/24) TC  snX4*nyROX*( I*9Y8PX*DYPDX/3 )/24 
CYTFRMs-0.500*rS0Y*ALMY™*0.  500*'"  SO  Y2*PXBOY-C  SO  Y3*  ( 0 Y^DX+DI  AG/  ( 1 2* 
10X1  ♦ 9LNYTM/24)  ♦CFnY4*OX«nYT(UDX?DY*PXPDY/3)/2't 

rALC  INDIV  SEPIFF  EPF  SFLF-INTFP 

CINTEPC1  tl)*-2/CS*<CXTFRM*rYTF5v) 

C0MPUTF  S E L F - 1 Mr  F 0 A CT  IrN 


no  220  I S*l*  2 

XS*(FLOATIIE|-l .5D9)*nx 

do  ?po  js»i,2 


LOOP  TO  CALC  ADJACENT  p/tTCH  TN^EP 


IF  I IC^JS.EQ.U  GO  TP  22' 
Y$MFLf'AT(  JS)-1 .500)*OY 
0?  210  I TNT*1 , INTPTC 
XP=FLn*T ( 1 1 NT-1 ) *DX  TNT 


NUMFF  TNT  WRT  X irnP 


X2TM2=XS*XP 
X2™2*X2tM2*X2TM2 
O',  ? 00  JIN-sl  f TNTPTS 
YP*FJ_nAT  ( JJNT  -1  ) * 0 Y I T 


Ml i m f p INT  Y Lrno 


Y?TMr  YS  *YP 

0 *0S  OPT ( X2Tm2+Y2TM*Y2TH  1 
“ I NTY(  jt  Nt)  *cof xp ( -?*C S*R ) /p 


EVAL  IKTFr,F*N,3 


200  CONTINUE 

'.ALL  0WE90LC  INTY,INTPTStrYINT,riNTXUTNT)) 

C INT  WR  T Y Tr  YIEL"  X T NTFfjo  a NO 

210  CONTI NIJP 

CALL  OWFOOL("TNTX,TN-rPTS,nxTNT,CIMTFP(TS,  JSI  » 

C.  INT  WPT  X 

220  CONTTNUF 
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02050 

02060 

02070 

02090 
02000 
02100 
02110 
02120 
02130 
02140 
02150 
02160 
021  70 
021  80 
02190 
022nn 
0?’10 
02220 
02230 
02240 
02250 
02260 
02270 
o 22  8" 

02290 
02  300 
"*2310 
0 2‘,?0 
02330 
0234' 

023  50 
02360 
02370 
02  3 80 
02390 
02400 
02410 
02420 
02430 
02440 
02450 
0 2460 

024  70 
•02460 
02490 
02  5 CO 
0251C 
02520 
02530 
02540 
02550 
0256'’ 
02570 
02580 
02590 
92601 
0261  C 
"*2620 
02630 
02640 
02650 


UlkUmttAi 


in  r>  >o  i i r»  .1  ri  n i n “>  ct  r>  -j  i ^ r»  r»  ->  ->  i ,t  n 


0"  ?50  T«*1fMXT2 
X2TM?«r>FLr'AT(  fC-i  ) #f)x 

O''  250  J5«1,MYT2 

Lonpc  FnP  FFMMMPPR.  P F INTFPACHO 
M CALCF"  BY  rNF-oT  pFFTAN0  PIILF 
!P| IS«-JS.LT .^.PR.TS.Fn.?.AMO.J$ . FO ..  2 ) PC  TP  250 
Y2TM=F|.nATt  JS-U*0Y 
3=0S0PTl X2TM2tY2tm*Y,tv) 

CINTFP(TS,JS )=PDFXP(-2*CP*P  )/o*OX*DY 

250  CONTINUE 

FNO  OP  LOOP  TP  FILL  I toTFP  * P T I Pi  MMPIX 

5PGTM  PnA'STPUPT!PN  nc  MPMFNT  MATpix 


O'1  350  I M=1  t MX 
03  ’50  JM*l *NY 


I*(  JM-n*NXMM 


MSU  »S  = 0 

Dn  330  JS*l.NYMl 
jni * TABS  f J«-JS )«•! 


J02=JMfJS 


OP  3\  0 I S=1  , MXM1 

toi= T a as ( is-Im  )*i 

TO?rI«;+TM 

J=( JS-l )*MX*TS 


T M D F X T N 0 PF  Mi Tf  H-onI NTS  OVFR  EMT 
T P F OlJArc«NT 

PNF-PtP  Mi.TfM-PT  tnofa  gfm*fd 
CPL'WTCF  ALONG  x-niPFCTlnN 


TNPFX  PV FF  SPUPPP  0 at PHFS  Y-OI® 

1ST  AMP  2MD  Oil  AD  J 'PIFFEPFMPF 
INDFX' 

?on  r,  4TM  Qi;ip  J 'DIFFFPFMrF 
TMOFX' 

v|ptp  ^hAt  'PIFFFPFN'F  INDICES'  »P 
F = i tnpfx  01 FFFp  FNr F ' H FOP  TPF 
SAKF  PF  P n P T F i v)  INDEXING 

INPFX  PVFF  Sr-UBCF  ''»t,‘hfc  X-OIF 

1ST  r,  4TH  OIIFO  1 i J FF  • OP  X 1 

2MP  F.  xrrj  0|i?0  • riTPr  rfjOF-X' 


OVF-PJM  snjQFF-PT  TMOFX 
COsC  I NT  FP  ( I 01  , JD1  ) ♦NSM  TT  T*r  TKrpp  ( ID?, JD?) 

SI iv  pF  jnijprp  pp NT  FROM  01  f.  OH  I 
CFsNS'HT  ♦CI'ITFP  (t  O’, JO’  )-+M  SM  T V*C  I NTFp  ( ! or  , J0  2) 

fjm  cf  source  :"n,t  pccm  on  r.  orv 

:''i’xH,j-NPij«s)=ro«-rF 

SUPMAT  FNTRY  Fop  X-OT  p CijrF'S 

C^i ty ( I * J ) =cn-r  f 

$HRM/iT  EMTCY  FOP  Y-Cto  CUPP'S 
NOTE  THAT  F Vc  M O'  s FONT  MFP.ATTVF 
pop  Y-OIP  CUPP'S 

310  OONTTYUF 

FMO  ’F  SnUpPF  LpnP  FOR  imt^ptop  PASHES 
• Z 0NC  T PUr T T OM  OF  snuorc  ttomc  FC  C M S’C(X)«i  EDGF  FPL  LnW  S 


TD1=T'BS(NX-TM)*T 
T02  = NXT-IM 
J* J S*NX 


. -v:  U/-X..  V.  ■ J.1*. 


n«rtMT«(  TOT  « JO\  T*rTNTFF<  jnpt  j(?  2»  03?  70 

r SUM  PF  «?UPTE  CPN^  FC0M  (Jl  6 QIII  0328C 

IK  «NS^T  T*2T  MTFB  ( I O’  , J01  )4.f'SMTV*riNTFPn01,  JD?)  03?90 

c stt*  pf  sruDfF  cpnt  ffom  on  r.  oiv  03300 

P "V( I f J J * ( T P-CF )*FOGPAP  03310 

C S'JBMAT  FMr  Y FPP  Y-OTP  rlJRF'S  03320 

r smc  that  fvn  0'S  PPM  NFPAT’VP  0333n 

* FPR  Y-PTP  PIJPP  • S 03340 

I'M  JS.NF.JPGF«smt»AF«-n  > p,p  TP  325  03350 

NF'JF'-  MINCMNSUPS  + I.NOBFJmI  03360 

3P  tp  330  03370 

315  PM*TX{  T,  J-MSUAS  )-'PF*rr  J *Fnr,?«c  03380 

c SJPMAT  FMTRY  Fn»  X-QIP  CUSP'S  03390 

p 03400 

P PiJO  PPUTIMF  FPC  ABS(X)=*  PPGF  tfgms  03410 

P 03420 

330  rTNTTM'JF  03430 

P 03440 

C PNO  LPPP  PVF°  Y PPOPP  FPR  IMTf: jpo  PfTHFS  P345C 

P 03460 

f PPGTN  ROJTTVF  Fpp  Pp»cTOllr  tt  PS)  fF  FPIIFP  F t^oms  fpp  ABMY)=B  EDGF  03470 

C 03480 

J01=T«BS(NY-JM)+1  0349C 

JD2*MY*JM  03500 

NSUPSJ=MR(JPS  03510 

V>IJBS=0  03520 

CP  340  T$«l » N!  X**l  03530 

C IMPFX  PPWN  X »"P?D'P  1MFRTPR  03540 

r °t tp HF  p O’ 5 50 

TP1=  T 3 3C ( TS- I M ) +1  03560 

T02=IS+TM  03570 

J = ( NYMl ) *NX*T S O’ 5 80 

r?sriNTFR(  101 , JD1  I OKMJT  I*PtMTFf  ( ID2,  JD2)  P359 r 

1 C'F  S rURr F PPNT  FPPm  oi  c 01 II  03600 

:P  = HS«TT  *FTNTFR  ( 105,  J01>  ♦NFMTV*PT^TEFUP1,  JO?)  03610 

c SUM  n f c nlJR  P F FPf'T  FPPM  on  r.  PIV  036?P 


335 


r 

c 


340 


r 

r 

r 

c 

r 


C*ATX(  I,  J-MSUBSJ)  =(PF  + rT')  *POP,F/'P 

SUP  MAT  F-'TPY  FPP  X-OJP  CUFP'S 
TF|  TS.MP.TOPF4S  fA'Sll°FTl  n F’P  TP  336 
\JSUPR  = MINO(  MSUnS+l,NPRFI  M ) 

OP  t-i  340 

* M*  TY  ( I , J-MFUF.F)  = (PP-PF)  *FOGFAf 

SUPM^T  FMTP Y FPP  Y-OTC  PiJPR'S 
sj-'Tp  thAt  PVFH  O'F  PC  NT  NEGATIVE 
FPP  Y-PIP  PlIPP'S 

PPNTINiJP 


FNO  ® 'UT1  Nr  PPP  ABS(Y)=B  Fpr.r 
P T*M  F TP  IIP  T I PN  PF  F nRN  co  PPJPCF  TEFM 


TP1=T4BS(NX-Im )♦! 

I0?=K'X4-TM 

Jr NX*MY 

rp=riNTcoi  ini.joi  )+ns«tt  ro'TVTFcn  p?f  jc?) 

SUM  PF  SiPURFF  C.PM  FPPM  01  t OUT 
CF*NSMn*-lM’-eP(lD2,J01  ) ♦NSMlV*rTNTEP(IPl  ,JO?) 

RUM  C’F  RPUBCE  PPM  FPPM  on  f.  OIV 
TF(NY.NP.JPPE7S(NP°cJ))  PMfTX(i , J-NPBFjM)=(rF+fP)*FDG5 

SUPMC'  bntry  ppp  X-DCP  CUPP'S 
Ic(  N'X.NE.TPFE  AS  (NPPFt  ) ) P’^-rY(  I,  J-NPCFTM  ) = I FC-P  F ) * FOG? 


03630 

03640 

03650 

0366C 

03670 

0368C 

03690 

03700 

03710 

03720 

03730 

03740 

03750 

03760 

03770 

03780 

03790 

03800 

03810 

03820 

038’C 

03B40 

03850 

9386r 

03870 
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JL 


r 

n 

r 


r 

r 


Stinv.si  Ff.'TP  v Ft  ? Y-pT  > riic'jic 
\'37F  "MfT  rVFf  ?*c  rr‘~  \Fr;/TJvF 
FP  R Y-^J0  fU rP'S 

350  rn^TTMUc 

FSO  OF  MPM  = KJT  M»TOTX  T^TPFAFTJriV  r P* ST°  Uf  'r  I CH 

pfo  r n poiit!mf  tp  fntfo  HPMPr-FroFnjs  FPL *n  fxp  ans  icn  rni  'S  in  m*trix 
360  Npes=2*NPPE 

HZGHFSt  prdff  bfc,fcL  FUNPt  TPN  IN 
HPMrrFKJFPDS  RPL'N  EXPANSION 


r 

r 


CNF  LESS  IF  EVEN  INDFX  EXPANS ICM 


INDFXUG  TH»II  vATH-t>TS 
PnL/*P  CPPPP'S  °F  MATPH-PTS 


364 


368 


IF(N!NDX.C0.2  » NRF<=NPCR-1 

O-i  400  IM  = 1,NX 
X=(FLnAT(TM|-o.5np)*nx 
Op  400  JM  = 1 « NY 
Y=(FLPAT ( JMJ-0.S00)*PY 
I=(  JM-1.  )*NXTIM 

PH1=98T4N(Y/X ) 

RHO=OFOPT(X*X4-Y*Y) 

OR8Fr.  = ?*PIM8r,  (CS  )*PHP 
f)TM89  0s-2*0BC6L(r:S>*eHp 

APOU^FN’-  of  PESSFI  PN'I 
1 c ■.  P AB5  C 01*  f PG/  On  Ap  G)  . t * ♦ 1 . F-2*  I GOTO  364 

IF  r Ft  I 4p0  SKIP  TO  F FAL  EES  CALL 
PK  LI  8S0  J Z(  no  Ap0,  OT  MAP  r,,  Or'RC<:.  n^PFS.  VP  FS.C.  PO,  16,  IP  PP  ,rUMl,D'JM2,D 
1U*3  .0IIM4) 

G'T  TAPIF  pf  PESSEL  FlIIPTIPMS 

r,p  T0  368 

-lU  PSLJZIOR  APG,OPncr  ,M°ce,  0.00,16,  I FTP  , oijm  1 , 0UM?  ) 

<-f  Lt  ZPPOZ(  01  mpFS  ,,•>*(  NPf«4-’  ) » 

POLL  ZEPpZ(OTMPFP , 2*(N0ESt1 M 

SFT  |ir>  PUFF  FOJL  P F S FI  IMFT  IrN  S 

rritry  |i  |sn 
r'IVTM(l  1=0 

ZFRp  !ST  t£cm  rpFF  p PI STRUCT TPN 
VEFTPP  c 

I VOE  X (J  rzi  r nF  rrPF  r nNF  t p 
V^TT'-F 


0P  37p  1 1 =1  ,MPFFP1 

TM0X  = 2*n-MTA'OX 

1FJ INOX.FO.P)  or  to  370 

AC  G=nFLOA'r(  I VOX -l  )*PHT 

.APODMFIT  pf  sin  fn 

rpcs  = nrMPLX(  DPP  f«  ( T*'PX  ) , PTMPFS(  IN'DX  ) ) 
rrncrv(Ti)  = proF(APr,  )*'PFS*4<'r'I 
*STN™(  II  | = D S T N ( AB0)*rncST4*PI 


'•Air  SPFIFS  INDEX 

sktp  rfir  of  Prirw  tfcn  Fpp  zero 
I N0F  X - IT  WAS  SET  Tr  ZEP°  ARPVE 


370 


- OMTT  M|JC 

rp  280  JJ  = 1,N|PPFJ 


J= JPPPf S(JJI*NX 
T NO  X=  ?* J J-NT NPX 


rtir  rrpFF  C p N 5 T P Ur  T T r K TFPMS 

(_-r-p  Tf.  c F P L A 0 E PPL  ’ S FPP  PF  E * S SI 
Gk:  ED  J T FP  MS 

IVPFX  P F COL  PF I NO  PPPLACFO 
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02980 

P3BOC 
03400 
03913 
03920 
039’n 
93940 
03950 
0 39  63 
0397r 
03980 
03990 
94000 
04310 
04323 

040  30 
04040 
049RP 
04360 
04970 
04,9  80 
04090 
04109 
941  K 

041  20 
04133 
04140 
34150 
0416r 
04170 
0418P 

041  9p 
042C0 
04213 
94?2p 
04230 
9<*243 

042  50 
04269 
042  70 
04?  80 
04290 
04303 
04310 
04320 
04330 
04340 
04350 
p43  6p 
04370 
04380 

nt,iqn 

04400 

04410 

04420 

0443C 

04449 

04450 
04463 
94473 
0448C 


Btffttfilliiizriiii 
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I 

I 

i 


; S'ctp?  I'vrtX  ftp  prpi^rj^r,  t^c-m  04499 

0?  t-i  ( 371  *37?)  ,NSr  ("*  94500 

SFl Fr  T Pacpf1?  S F P T E c fr  F r errnpoi  9451C 

SO  Tf'  Y 0 YM mp Tp y r 0'in  I T I ~fJ  9452'* 

NS(*T=?  I NCi  I r .*T  FS  rrc  IKFS  IM  X 94590 

0  HF  r F * T ro  94540 

371  rHnv(i,jj)  = -Pi/(2*rsi*(9.or, i .pr  >**  t*'ox*i  r«i  ntm  { j j j _r  c tf  tm  < j J+ u ( ^4550 

rHiM(  vjorr  t , JJ ) =- DJ  / ( 2 *r  c ) *(9.  Q9  ,1  . HO  )**  \ N|OX  * lOOP  ( JJ  1 ♦ 94560 

lCr95TM( JJ+1 ))  94570 

r,n  m 3pn  94580 

372  "|-nM(ttjj)=-D!/(j*r?)*(O.POtl.P0)*MH0X+1  >*trrr««TM(  JJ-*-’  )-  945  99 

1 rr"*f  tm(  j j j 1 ''4690 

■H9M(T*NPT'-H<;  ,JJ)  = -pt/  <?*<-«  1*10.09,  ’ .OP  )**<  IF  DX*1  1*TC  INTM1  JJ  )♦  94610 

KSTXTNM  J J + \n  946  29 

380  C°NTIMl)P  ''463'’ 

On  390  I T =1  ,MPPFI  94640 

J=(MY-1  )*MX*!F>PF»S(TT  l^'PT-wc  94659 

LTD  T r ’FPl.TF  fPL'*  FT  PCFAF;,)  ''466n 
J J = T T ♦NPRC J 04670 

0‘JPD  T TF-M5  94680 

T M 0 X = ? * ( I T+NPPEJ)  -VTNDX  ''46  90 

09  T9  < 381,38’) ,M$r T 0470C 

381  OH?M(  l,jj  |s-P!/(?*o?)*(9.no,  i.DO)**!tf:nx*{f:STF,TM(TTfNDPFJ)-  94  710 

10OTF't^(TT*mppcJ  + 1 ))  04720 

0M9M1  I*MPt:hS  , JJ>  =-D?  / 1 ’ *r  5 ) *10. 00 ,1 . 001**1  \DX*lrFrST>'(  II  +F'PF  C J)  ♦ 04730 

lCr9S7M( II+MPOFJ*l  1 1 94749 

30  T9  390  04750 

38?  0HOM(  j tjj)  = -pT/<;?*rc)*(c>.H|o,  1.301**1  IMnx+lJMnCOSTMlTT^  PkFJ+l  1-  "<47  69 

irrnp*M(IT*\'POFJ)l  947  7r 

r.H’,M(T*M0TrH5f  jjj-.or/i ?*rc)*( Q.oo, l.P9j**(  imox  + 1 J*  94780 

1  ( 0 5 I M™  1 T t*\|PPFJ  » *0  c t m TM 1 T T fMPorjt  1 J ) 947  9.9 

390  CONTINUF  ''480° 

400  “ONTTNUP  04310 

r 04329 

r rtgr)  nF  MCMFF’T  M/IJPTX  fnv  c T n t IF”  I rt.'  "4031* 

r 04340 

405  onUTIMlJP  043  59 

0«LI  5PRH0M(CM*TX  ,NPTOMF  r TMF -F'PR  F j , r nTMif  ,inj  Ml , -40 a<- 

1 TMfTY.Npri-MF  , MPT'-HC -NPP  F I ,F  ni  Ml  , \L  TMl  , 04J7r 

2 rHrH,Mn»Mry,Mr1TMr  j,rM«T<,,^  ri  MF  J , r.  n j mj  ,r  pFT  j 04800 

Fo»-=rnsBS(rM«TX(  t tl  j j o48or 

noLATC  = Or)PT/PPAT  94900 

9PI-F16,  29)  r5HNn3,rPl.ATF  04910 

20  PnOMATJ*  • ,5X,2P12.4,5X,2F12.41  04920 

RPTURU  94939 

FM9  94640 
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SUBROUTINE  $PPH0MICMATltNIl,NJl,NDIMlIfNDIMlJ,CMAT2tNI2tNJ2,NDIM2I  00010 
1 » NDI  H2J  *CMAT3f  NDI  M3 I » NDI  M3  J t CM  AT4 ; NDI  M4 1 » NDIM4  Jf  CDET  ) 0C020 

IMPLICIT  COMPLEX*  16(C)»REAL*8(  A »B» D-H  »C-2 I 00030 

00040 

SUBROUTINE  TO  OIAGCMLIZE  ANC  CALC  DETERMINANT  OF  A SPARCELY-  00050 

COUPLEC  MATRIX  00060 

BY  L H PEARSON  7/74  CC070 

RE VI SEO  5/75  00075 

00080 

01  PENSION  CM ATI  (NC IMl I , N DI Ml J 1 # CMAT2 INDIM2I. NDIM2J I • CMAT31 NDI M3I , CC090 

1NDIM3J) »CNAT4(NDI M4I»ND! M4 J)  00100 

NI3*NIltNI2  00110 

NJ3*NI3-NJ2-NJ1  0C12O 

CALL  ZEROZ (CMAT4f  4*N0I M4 I*NOI M4  Jl  00130 

CDET  «l  0C140 

INITIALIZE  PRODUCT  ACCUMULATOR  CC150 

NPR-3 

N J 1 M 1»  NJ  1- 1 CGI  60 

NJ1L-NJ1 

IF( NJ2+NJ3.GE . 1 ) GO  TO  95 
NJ1L«NJ1L-1 
NPR*  l 

95  DO  155  M* If  NJ  1L  00170 

INDEX  ACROSS  COL  001 80 

HP1-M  + 1 00190 

FMAX-CDABS<CMATl<MtMI)  CC2C0 

K»M  00210 

IFJMP1.GT.NI1  ) GO  TO  105  00220 

DO  100  I *MP1 (Nil  00230 

LOCF  TC  SEARCH  FOR  PIVOT  IN  MTH  00240 

COL  00250 

FCK*CDABS (CMATllItMII  0C260 

IF  I FCK.LE.FMAX)  GC  TO  100  00270 

K=  I 00280 

C IF  LARGER  ELEMENT  FOUND  MARK  ROW  00290 

FMA  X*FCK  00300 

C USE  NEW  LARGE  ELEMENT  AS  CCMPARI-  00310 

C SON  VALUE  C0320 

100  CONTINUE  00330 

105  CSTOR*CMATl( K»M ) 00340 

C SAVE  VAL  OF  PIVOT  ELEMENT  00350 

CDET»CDET*CSTOR  00360 

C MULT  PIVOT  INTO  PROD  ACCUMULATOR  CC370 

IF (K.EQ. M)  GO  TC  115  C0380 

C IF  PIVOT  CN  DI  AG  SKIP  ROW  EXCF  00390 

COET  =-CDET  00400 

C CHANGE  SIGN  BECAUSE  OF  ROW  EXCH  00410 

107  DO  110  J«Mf NJ  1 00420 

C LOOP  TO  EXCH  DIAG  AND  PIVOT  ROWS  00430 

CSTC*CMAT1 (Kf J)  00440 

CMATK  K,  JI«CMATKM,J)  00450 

CMAT1 ( M»  J) *CSTC  CC460 

110  CONTINUE  00470 

IFINJ3.LT. II  GO  TO  115  00475 

DO  112  J>1»NJ3  CG480 

CST0-CMAT3I K,J)  00490 

CMAT3IKfJ)*CMAT3(MtJ)  CC5C0 

CMAT3( Mf JI-CSTC  00510 

112  CONTINUE  0C520 

115  CONTINUE  00560 

IF ( MP 1 .GT.NI 1 1 GO  TO  155  00570 


1 


. 

\ 


8 

i 


53 


MBip1 


jlPfl  UJM| 


T 


i- 


00  150  I-«P1»M1 

00560 

1 j 

c 

ELIMINATION  LUO''  FOR  CHAT  1 

00590 

CFAC*CMAT1< It HI/CSTOR 

0C6C0 

1 

c 

EL IM1NATI0N  FACTOR 

00610 

IF  ( HP  1 .GT.NJ 1 ) GO  TO  125 

00620 

00  120  J"HPlt  NJ1 

00630 

3 

c 

LOOP  ACROSS  ROM  IN  CHAT3 

00640 

'j 

CMATHI,  J)«CMAT1(  I,J)-CMAT1<  M,J)*CFAC 

00650 

■ \ 

120 

CONTINUE 

00660 

IFINJ3.LT.il  GO  TC  150 

00665 

• | 

12S 

00  130  J*ltNJ 3 

00670 

C 

LOCP  ACROSS  ROM  IN  CHAT 3 

00680 

1 

CMAT3I  It  J)-CMAT3(I  , J)-CMAT3C  Ht  J)*CFAC 

00690 

; 

130 

CONTINUE 

00700 

150 

CONTINUE 

00720 

155 

CONTINUE 

00730 

V 

M4-NI1-NJ1 

00740 

] 

IF ( M4.LE.0I  GO  TC  290 

00750 

C 

00760 

> 

1 

C 

00770 

H 

C 

BEGIN  ROUTINE  TO  CREATE/ 'DIAGONAL! ZE  • CHAT4 

00780 

f 

C 

00790 

si 

NPIV-NI4 

008C0 

‘i 

IFINI4.GT.NJ2)  NPIV-NJ2 

00810 

DO  250  M-l.NPIV 

00820 

i 

C 

INDEX  ACROSS  COL  FCR  CM AT 4 Cl AG 

00830 

V 

MP  l-M*l 

00840 

FHAX-CCABS(CMAT2(1,M)) 

00850 

i 

K«*  1 

00860 

V 

• i 

IF(NI2.LT.2I  GO  TO  205 

00870 

■f 

a 

00  200  I>2tNI2 

00880 

C 

LOCP  TC  SEARCF  FOR  PIVOT  IN  MTH 

00890 

C 

COL 

00900 

0 

FCK«C0ABSlCHAT21ItM) 

00910 

J 

IF(  FCK .LE.FMAX)  GC  TC  200 

00920 

s 

K«I 

CC930 

C 

IF  LAPGER  ELEMENT  FOUND  MARK  ROM 

00940 

■ 

FMAX-FCK 

00950 

•i 

C 

USE  NEW  LARGE  ELEMENT  AS  COMPARI- 

00960 

r 

C 

SON  VALUE 

00970 

200 

CONTINUE 

00980 

'■ 

205 

CST0R*CHAT2IK«M) 

C0990 

* 

!> 

C 

SAVE  VAL  CF  PIVOT  ELEMENT 

01000 

,* 

k 

COET-CDET*CSTOR 

01010 

'* 

C 

MULT  PIVOT  INTO  PROD  ACCUM 

01020 

CDET—COET 

01030 

C 

CHANGE  SIGN  OF  DETERM  BECAUSE  OF 

01040 

4 

C 

EXCHANGE  FROM  CMAT2  TO  CMAT4 

01050 

i 

DO  210  J»M,NJ2 

01060 

jjfl 

C 

LOOP  TO  EXCHANGE  01  AG  AND  PIVOT  R 

01070 

1 

C 

ROWS 

01080 

1 

CST0*CMAT4(M, J) 

01090 

■m 

CMAT4C M, Jl»CMAT2(KtJ) 

01100 

I 

CHAT2(K, JI-CSTC 

OHIO 

210 

CONTINUE 

01120 

1 

K3-K4NI1 

01130 

| 

M3-NJI+M 

01140 

S 

IFINJ3.LT. 1)  GO  TO  213 

01145 

'j 

00  212  J«ltNJ3 

01150 

m 

CST0-CMAT3(K3tJ  ) 

01160 

1 

CMAT3(K3tJ)«CMAT3(M3tJ)  54 

Oil  70 

1 

O <1  O O O O o o 


212 

213 

235 


C 

C 


240 
2 42 


245 

250 


w 


CMAT3CM3,  JI-CSTC 
CONT  I NOE 

IFIM2.LT. 1)  GC  TC  250 
DO  250  I ■ 1, NI 2 

LOOP  TO  CARRY  ELIMINATION  INTO 
CMAT2 

I3*N 11*1 

CFAC«CMAT2II,M1/CST0R 
IFC  HP  1 .GT.NJ2 1 GC  TC  242 
00  240  J*MP1,  NJ2 

LOCP  ACROSS  ROM  OF  CM AT 2 
CMAT2( I, JI-CMAT2CI ,J)-CMAT4CM, J) >CFAC 
CONTINUE 

IFCNJ3.LT.il  GO  TO  250 
00  245  J*1 » NJ  3 

LOOP  ACROSS  ROM  OF  CM  AT 3 
CM  AT  3 ( I3,J)»CMAT3< 13, J I-CMAT3C M3, J)*CFAC 
CONTINUE 
CONTINUE 


END  ROUTINE  TC  ' Cl ACON AL IZ E • CMAT4 


290  IFCNI4.GE.NJ2 1 GO  TO  350 


IF  DIAGONAL  DOES  NOT  PASS  THRU 
SKIP  CIAGGNALI2ATICN  FOR  CM  AT  2 


BEGIN  ROUTINE  TO  'DIAGONALIZE*  CMAT2 


t 

3l 


t 


295 


C 

C 


C 


C 

c 

300 

305 


C 


C 

C 

C 


NI 4P 1=N 14*1 
NJ2L-NJ2 

IFCNJ3.GE.1 1 GC  TC  295 

NJ2L*NJ2L-1 

NPR*2 

DO  350  M*NI 4P  1 , N J 2 L 

MI*M-NI4 

M3*MI*NI  l 

MP1=M*1 

MIPl»MIU 

FMAX*C0ABSCCMAT2CKI,M)  I 
K«MI 

IFCMIP1.GT.NI2I  GO  TO  305 
00  300  I *MI PI  ,N 12 


FCK*CDAB5CCMAT2CI,M) I 
IFCFCK .LE.FMAXI  GO  TO  300 
K*I 

FMAX«FCK 


CONTINUE 

CST0R«CMAT2CK,M) 

K3«K4NI1 

CDET«COET*CSTOR 

IFCK.EQ.MI I GO  TC  315 

COET—CCET 


LOCP  TO  SEARCH  FOR  PIVCT  IN  MTH 
COL 

IF  LAFGEP  ELEMENT  FOUND  MAPK  ROW 

USE  NEW  LAPGE  ELEMENT  AS  COMPARI- 
SON VALUE 

SAVE  VAL  CF  PIVOT  ELEMENT 

MULT  PIVOT  INTO  PRCO  ACCUMULATOR 
IF  PIVOT  ON  DIAG  SKIP  POM  EXCH 
CHANGE  SIGN  BECAUSE  OF  ROW  EXCH 


01180 
01190 
01225 
01230 
01240 
01250 
01260 
C1270 
01280 
01290 
01300 
01310 
C1320 
01325 
C1330 
01340 
01350 
01360 
01380 
01390 
C1400 
01410 
01420 
C 1430 
01440 
01450 
C1460 
01470 
01480 
01482 
01484 
01486 
01488 
01492 
C1500 
C 15 10 
01515 
01520 
C 1530 
01540 
C1550 
01560 
01570 
01580 
01590 
01600 
01610 
01620 
C1630 
01640 
01650 
01660 
01670 
01680 
01690 
01700 
01710 
01720 
01730 
01740 
01750 


\ 

5 


.d 


"Tirpn)^ 


310 


c 

c 

c 


c 

C 


C 

C 


400 

405 


C 

C 

C 

C 


LOOP  FOR  EXCH  IN  CMAT3 


312 

315 


ELI  PI  NATI  CN  LCOP 


330 

335 


345 

350 


390 


00  310  J-M,  NJ2 
CST0-CMAT2CK, Jl 
CN AT2( K«  Jl«CNAT2t  Ml*  J1 
CMAT24HI, Jl-CSTC 
CONTINUE 

IFCNJ3.LT. 1)  GO  TC  315 
00  312  J«1»NJ3 

CST0*CMAT3CK3  * J) 

CHAT  3( K3fJ5»CMAT3CM3tJ) 

CMAT3CM3,  Jl-CSTO 
CONTINUE 
CONTINUE 

IF ( MIP1 .GT .NI2 ) GO  TO  390 
DO  350  I-MIP1.NI2 

I3-I4NU 
CFAC-CMAT2C I , Ml /C STOR 
IFCMP1.GT.NJ2)  GO  TO  335 
DO  330  J«MP1,NJ2 

LCCF  ACROSS  ROW  OF  CMAT2 
CM  AT  2(  I,  J1-CMAT2C  If  J1-CMAT2CMI,  J)*CFAC 
CONTINUE 

IFINJ3.LT.  11  GO  TC  350 
DO  345  J*l , NJ3 

LOOP  ACROSS  ROW  IN  CM  AT  3 
CMAT3C I3,J)*CMAT3( I3f J1-CMAT3I  M3,J)*CFAC 
CONTINUE 
CONTINUE 

BEGIN  ROUTINE  TO  * Cl  AGCNAL I Z E • CMAT3 
NJ3M1-NJ3-1 

IFINJ3Ml.LT. I 1 GC  TC  455 
DO  450  M*l,NJ3Ml 


INDEX  ACROSS  COL 


MPl»M*l 

MI  *M>NJ 1+N  J2 

MIP1-MI+1 

FMAX»C0ABSCCMAT3CRIfM> ) 
K«MI 

IFCMIPl.GT.NI3)  GC  TO  405 
DO  400  I*MIPl,NI3 


FCK-C0ABSCCMAT3C I ,M) 1 

I F ( FCK  .LE.FMAX  1 HO  TO  400 

K*I 

FMAX«FCK 

CONTINUE 

C$TCR*CMAT3C  K,M) 

CDET  *COET*CSTOP 
IFCK.EC.MI)  GC  TO  415 
C0ET»-C0ET 


LOCP  TC  SEARCH  FOP  PIVCT  IN  MTH 
CCL 


IF  LARGER  ELEMENT  FOUNC  MARK  RCW 

USE  NEW  LARGE  ELEMENT  AS  COMPARI- 
SO  VALUE 


SAVE  VAL  CF  PIVOT  ELEMENT 
MULT  PIVOT  INTO  PPCC  ACCUMULATOR 
IF  PIVCT  CN  DIAG  SKIP  ROW  EXCH 
CHANGE  SIGN  BECAUSE  OF  ROW  EXCH 


Cl  760 
01770 
01780 
01790 
01800 
C1805 
01810 
01820 
01830 
01840 
C1850 
01860 
01900 
C1910 
01920 
01930 
C 1940 
01950 
01960 
C197C 
01980 
C199C 
02000 
02005 
C2010 
02020 
02030 
02040 
C2060 
02070 
C2080 
02090 
02100 
02110 
02120 
02130 
02140 
02150 
C2160 
02170 
02180 
C2190 
C22C0 
02210 
C2220 
02230 
02 240 
02250 
02260 
C22  70 
02280 
02290 
C2300 
02310 
02320 
C2330 
02340 
C2350 
02360 
02370 
C2380 


56 


I 

'• 


I 


00  410  J-MVNJ3 

02390 

c 

LOOP  TO  EACH  OIAG 

ANO  PIVCT  RCWS 

02400 

CSTO-CMAT3IK. J| 

02410 

CMAT3(K# J)»CMAT3(M1,J) 

02420 

CMAT3(HI,J)-CST0 

C2430 

410 

CONTINUE 

02440 

415 

CONTINUE 

02480 

00  450  I-MIP1,NI3 

02490 

C 

ELI MINATI CN  LCOP 

02500 

CFAC*CHAT3(  ItM)/CST0R 

C2510 

00  445  J«MPl»NJ3 

02520 

C 

LOOP  ACROSS  RCM  IN 

1 CMAT3 

02530 

CHMT2!  J,J)-CMAT3CItJl-CHAT3(Pl, 

J)*CFAC 

02540 

445 

CONTINUE 

02550 

450 

CONTINUE 

02570 

455 

GO  TO  (461 1462  »463 1 » NFR 

02572 

*61 

CDET>C0ET*CHAT1(NI 1.NJ1) 

02574 

RETURN 

C25  76 

462 

CDET*CDET*CMAT2  (NI2«NJ2) 

02578 

RETURN 

02582 

463 

CDET-CCET*CMAT3(NI3,NJ3» 

02584 

RETURN 

02600 

C 

MULT  LAST  ELEMENT 

I NTC  OETEPM 

02590 

END 

02610 

i 


1 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 

c 

81 


82 


83 

C 


C 

C 

C 

C 

C 

C 


C 

C 


SUBROUTINE  SRRSLVlCMATitN! 1. NJI, NDIM II .NDIMIJ, CMAT2, NI2 ,NJ2,NDIM2I 
l,NDIM2J,CMA?3,N0IM3I,N0IM3J,CMAT4,N0IM4l »NOI M4 J* CRHS,C SCLN) 

SUBROUTINE  TO  BACKSCLVE  A TRI ANGULAR  I ZEC  SYSTEM  OF  SPARCELY- 
COUPLEO  LINEAR  EQUATION 
BY  l N PEARSON  7/74 
REVISED  5/7S 


STORAGE  FORM  COMPATIBLE  WITH  THE  TRI ANGULAR I ZAT I ON  ROUTINE  SPARCE 

THE  ENTRY  'HCMSLV ' BELCW  ALLOWS  THE  SOLUTION  FOR  NATURAL  VECTORS 
OF  HOMOGENEOUS  SYSTEMS  PROVIDED  THE  DETERMINANT  OF  THE  SYSTEM  IS 
ZERO 


IMPLICIT  COMPLEX* 16(C) ,REAL*8( A,B,0-H,C-Z) 

DIMENSION  CMATlfNCIMl  I* NDIM1J ) , CHAT  2(NDIM2I«NDIM2J)  »CMAT3(  NDI M3I  ,N 
101 M3JI tCMAT4( KOI M4 1«NDIM4J )»CRHS(NOIM3I ) *CSCLN(N01M3I ) 

LOGICAL  LHOM 

A 

SETUP  FOR  I NHCMCGE  NECUS  SYSTEM 


LHCM-. FALSE. 

NI3-NI  1*NI  2 
NJ  3*NI 3-NJ1-NJ2 
NI4*NI l-NJl 

ND2-NJ2-M4 

NPR«3 

IF ( NJ3.LT.  1)  NPR-2 
IF ( NJ 3 *N J 2 .LT .1 ) NPR-l 


SET  INCICATCR  FOR  INHOM  ENTRY 
NO  ROWS  IN  COUPLING  SUBMATRIX 
NO  CF  COLS 

NO  ROWS  IK  SECONDARY  COUPLING 
SUBMATRIX 

NO  OF  DIAGONAL  TERMS  CF  MATRIX 
IN  CMAT2 


SET  INDICATOR  FOR  NULL  CMAT3 
DEGENERACY 

SET  INCICATQR  FOR  NULL  CMAT2  E 
CMAT3 


GO  TO  (81,82*83)  , NPR 

GO  MAKF  FIRST  DIVISION  FOR  RIGHT- 
MOST matrix 

CSCLN(NI3)*CRHS(NI3)/CMAT1(NI1,NJ1) 

GO  TO  IOC 

CSOLH(NI3)«CPH!i(NI2)/CMAT2(NI2,NJ2) 

GO  T 11  100 

CSOLNtNI  3)  *CR HS(NI3I/CMAT3(NI3*NJ3) 

SOLVE  FOR  'LAST*  UNKNOWN 


GO  TO  100 


GO  TO  SCLN  ROUTINES 


ENO  OF  SETUP  FOR  INFCM  SYSTEM 

BEGIN  ENTRY/SETUP  FOR  HOMOGENEOUS  SYSTEM 

ENTRY  HOMSL V( CMAT  1 , NI  1,NJ1 ,NDI Mil , NO I Ml J,CMAT2 ,NI 2 ,NJ2 ,NOI M2 1, NDIM 
L2J,CMAT3,N0IM3I , NO  I M3 J ,CM AT4 ,NC IM4 1 , NO  I M4J , CSOLN, NORO ) 

LHCM-.TPUE. 

LOGICAL  INDICATOR  FOR  HOMOGEN  SYS 


00010 
00020 
00030 
00040 
00050 
00052 
0C054 
00056 
00060 
OCOTO 
00080 
00090 
001 00 
00110 
00120 
00130 
00140 
00150 
00160 
00170 
00180 
00190 
00200 
00210 
OC220 
00230 
00240 
00250 
00260 
00270 
00280 
00290 
00300 


00320 

00330 

00340 

00350 

00360 

00370 

00380 

C0390 

00400 

00410 

00420 

0C430 

00440 
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nooo  or>  o <->  o o r>  o no  o o r»  o o o o o o 


NI3-NIHN12  00450 

NJ3-NI3-NJI-NJ2  00460 

Nl 4«NI i-NJ 1 00470 

ND2-NJ2-NI4  00480 

CS0LN1 NI 3)  *1  0C490 

ASSIGN  ARBITRARY  ELEMENT  IN  SCL'N  00500 

CC5I0 

END  SETUP  FCR  HCMCCCNECUS  ENTRY  0C52O 

00530 

BEGIN  BACKSOLVE  FOP  EQUATIONS  INVOLVING  ONLY  CHAT 3 (LAST  NJ3  EQS ) 0C540 

00550 

100  FMAX*CCA6S(CS0LN(NI3)I  00560 

IMAX*NI3  0C570 

IFINJ3.LT. 2)  GC  TC  2C0  00580 

SKIP  ROUTINE  IF  ONLY  LAST  VAPIABL  00590 

COUPLES  (IT  WAS  SOLVEO/ASSIGNEO  006C0 

ABOVE)  00610 

00  150  IC«2,NJ3  CC620 

ICMl*IC-l  OOl 30 

I«NI3~IC*1  00640 

I-M3-IC+1  00650 

C;L(  MATRIX  ROW  INCX  FROM  C0660 

COMPLEMENTARY  I NDX  00670 

J03*I-NJ1-NJ2  C0680 

COL  INOX  FCR  CM AT 3 WHICH  DEFINES  00690 
01  AG  CF  MATRIX  00700 

CSOM-O  00710 

00  110  J3C-1.ICM1  00720 

LOOP  TO  ACCLM  NEGATIVE  SUM  CF  00730 

PREVIOUSLY  CALC  1 D UNKNS  00740 

J3*NJ3*1-J3C  00750 

COL  OF  COEF  IN  CMAT3  00760 

J = NI3*1-J3C  00770 

ROW  OF  UNKN  IN  CSCLN  00780 

CSUM*CSUM-CMAT3 ( I , J 3 )*CSOL N( J I CC790 

110  CONTINUE  00800 

IF ( .NOT .LHOM ) CSUM*CSUM*CRHS(I ) 00810 

AOO  R H S TO  SUM  CCB20 

CSOLNl I)»CSUM/CMAT3(I , J03)  00830 

DIVIDE  BY  DIAG  COEF  0084C 

IF (CDABS(CSOLN( I ) ) .LE.FMAX  ) GOTO  150  C0850 

FMAX*CDABS(CSCLN( 111  00860 

IMAX* I 00870 

CHECK  FCR  MAX  ELEMENT  C0880 

150  CONTINUE  00890 

0C900 

BEGIN  ROUTINE  TC  SCLVE  FOR  ELEMENTS  INVOLVING  CM AT 3 C CMAT 2 00910 

00920 

20Q  IF ( NJ3 .GE.NI2 ) GO  TO  300  CC930 

C SKIP  ROUTINE  IF  DIAG  COES  NOT  00940 

C PASS  THRU  CMAT2  00950 

00  250  I C*  1 • NC2  CC960 

ICM 1*1  C- l 00970 

I2*NI2-NJ3U-IC  00980 

I3*NI3-NJ3U-IC  0C990 

JD2»NJ2*1-IC  01000 

NCMl*NJ3+IC-l  01010 

CSUM-0  01020 

IF(NJ3.LT.l)  GO  TO  215 

00  210  JC*ltNJ3  01030 

C LOOP  TO  SUM  CCNTRIB  FROM  CMAT3  01040 

59 


J3-NJ3+1-JC  01050 

J*»  N 1 3 ♦ 1— J C 01060 

CSUM»CSUM-CMAT3  I I3»J3MCSCLN(J)  01070 

210  CONTINUE  01080 

215  IFUCHl.LT.il  GO  TO  2?5  01090 

C SKIP  IF  NO  TERNS  CONTRIB  FR  CMAT2  CUOO 

DO  220  J2C-I.ICM1  OHIO 

J2*NJ2*l-J2C  01120 

J-M3-NJ3+1-J2C  Cl  130 

CSUM»CSUM-CMAT2U2tJ2l*CSCLN<JI  01140 

220  CONTINUE  01150 

225  IF ( • NOT.LHOM)  C SUM-CSUMaCRHS II 3 I CU60 

CS0LNII3)-CSUM/CMAT2(  12,  J02J  01170 

IF1C0ABS (CSOLNI  13  1 1 .LE.FHAX 1 GO  TO  250  01180 

FMAX«CDABS(C$CLN(I3I I 01190 

IMAX-I  ' 01200 

250  CONTINUE  01210 

C 01220 

C BEGIN  ROUTINE  TO  SOLVE  FOR  ELEMENTS  INVOLVING  CMAT3  6CMAT4  01230 

C 01240 

300  IF < NI4.LT. 1 1 GO  TC  400 

00  350  IC-I.NI4 

I4*NI4+1-IC  01260 

JD4*I4  01270 

I3«M  !♦  l—T C 01280 

CSUM*0  01290 

IFINJ3.LT.  11  C-0  TO  315 

DO  310  J3C*1.NJ3  C1300 

J3*NJ  3+ 1-J  3C  01310 

J-NI3+1-J3C  C1320 

CSUM»CSUM-CMAT3<  13,  J3 I *C SGuNI J I 01330 

310  CONTINUE  01340 

315  NSUBS-NC2+IC-1  01350 

C NC  CF  NCN-DIAG  CHAT 4 EL'S  IN  EO  01360 

IF( NSUBS .LT -1  I GO  TO  325  01370 

00  320  J4C*i  ,•  NSUBS  01380 

J4«NJ2*  1-J4C  01390 

J*NI3-NJ3U-44C  C14C0 

CSUM*CSUM-CMAT4(I4,J4I*CS0LN<J  1 C1410 

320  CONTINUE  01420 

325  IF( .NOT .LHCM1  CSUM=CSUM*CR FS II 3 1 C1430 

CS0LNU3)  *CSUN/CMT4<I4, 14  1 01440 

IF(COAeS(CSOLN(  13 1I.LE.FMAX1  GO  TO  350  01450 

FMAX-CCABS (CSCL  N 113  1 1 C1460 

IMA  X*I  3 01470 

350  CONTINUE  J1480 

C C1490 

C BEGIN  ROUTINE  TC  SCLVE  EC'S  INVOLVING  CHAT 3 C CHAT  1 C15C0 

C 01510 

400  IF1NJ1.LT. II  GC  TC  455 

00  450  I C*1 • N J1  01520 

I-NJl  + l-IC  01530 

ICMl-IC-1  C 1 540 

CSUM»0  01550 

1 r~  1 N J 3 «L  T . 1 1 GO  TO  415 

DO  410  J3C-1.NJ3  C 1560 

J3-NJ3M-J3C  01570 

J»NI3+1-J3C  C15e0 

CSUM-CSUM-CMAT3U  t J31*CSCLM  J)  01590 

410  CONTINUE  01600 

415  IFIICMl.lT.il  GC  TC  425  C1610 

60 


n n n M n ooor* 


DO  420  JC»1,  ICM1  £1620 

J-NJl+i-JC  01630 

CSUM-CSUM-CMATlUf  J)*CS0LNIJ)  01640 

420  CONTINUE  01650 

425  IF( .NOT.LHOM)  CSUM-CSU MKRFS (I)  01660 

CSOLNI I)*C  SUN/C NAT1I1 *1 ) 01670 

1F( C0A8S  (CSOLNt  I ) ).L  E»  FMAX  ) CO  TO  450  01680 

FMAX*COABS  ICSCLM  1 ) ) 01690 

IMAX-I  01700 

450  CONTINUE  01710 

01720 

01730 

ENO  OF  SOLUTION  01740 

01750 

455  IF( .NOT.LHOM)  RETURN  01760 

RETURN  IF  INHOM  SYSTEM  01770 

01780 

BEGIN  NORMALIZATION  ROUTINE  FOR  NATURAL  VECTCR  FOR  HGMCGENEOUS  01790 

CASE  01800 

01810 

CSCALE* l./CSOLNI I MAX)  01820 

00  500  I»l»  NI 3 Cl  830 

Cf>OLN(  I ) *CSCLM  I ) *CSCALE  C1840 

500  CONTINUE  01850 

RETURN  C1860 

END  01870 
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Tl~ r ■ Ciii iti  AinYfi  i?inr  •>w~ 


SUBPnjTTNg  ZFROZf IARRAY,N» 
3!HFNSnN  I ARRAY!  1 } 

Dn  100  T- 1, M 
I«RP«YU)»0 
CDMTfNUP 
RETyov 


09410 
09420 
094  3o 
094f0 
09450 
09463 
094  70 


i 


•! 


tl 

2 


I 


"B,^f 


""*■'  11 


100 


200 


SURRDUTTNE  OWEOOKFrN,  N,0FLTA,vlNT) 

IMPLICIT  REAL*8(A-H,n-Z) 

C0MPLEX*16  FCN,C,VINT 
DIMENSION  FCN(NI 
OIMFNSION  CORFU) 

DATA  COEF/2.DO»5.DO,l.nO,6.DO, 1.00, 5.00/ 

I F ( ( N-l I /6*6. EO.N-1 ) GH  Tr  LOO 
M*T  TF (6,11 

FORMATI'-OINCORPFrT  POINTS  n WEODLF*) 

• ■I/O 
CONTINUE 
VINT  *0 

DO  200  J*1 » N 
JC-»FF«J-((J-l  1/61*6 
VINT*VINT«-COEF(JCOFFJ*FCNI  J1 
CONTINUE 

VINT*  ( VTNT-FrN(  D-FCNCNI)  *(0.300, 0.00)*DfMPL  XI  DELTA  ,0.00) 
RETURN 


END 


09480 
09490 
095  OP 
09510 
09520 
09530 
09540 
09550 
''9560 
09570 
09580 
09590 

09600 

99610 

09620 

09630 

09640 

09650 

09660 
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... 


t n n n t t t 'i  n o o ,t  t n n ott  ii  m n i n o n n n T '(  o ^ n i n n o '.t  o n n lT  t t m o iT  n n o o ci 


I h 


C.Z  ANLVT 0, 

r. 

function 


USAGE 

PARAMETERS 


ZEROS  OF 
MULLER'S 


AN  ANALYTIC 
MFTHOO  WItH 


COMPLEX 


EPS 

NSIG 


- A 


- 1 


ZAN0967U 
ZAN0968D 
Z AN09690 
Z AN09700 
ZAN09710 
Z AN09720 
ZAN09730 
Z AN09740 
Z AN09750 
Z AN0976D 
ZANO9770 


KN 

NGUESS 

N 
X 


- A 


ITMAX 

INFFR 


DETERMINATION  CF 
FUNCTION  USING 
deflation 

CALL  ZANLYT  (F, FPS , NS IG,  KN ,NGUESS , N, X , ITMAX, 

INFER » T ER ) 

FUNCTION  SUBPROGRAM,  F(Z),  WPITtEN  py  THF 
USER  SPECIFYING  THF  EQUATION  WHOSF  ROOTS 
ARE  TO  RE  FOUND.  F MUST  RF  tYPF-NAMED  AS 
FOLLOWS  - COMPLEX  FUNCTION  F*16  IZ) 

ST  STOPPING  CRITERION.  A ROOT  Z IS  ACC EPT EDZ AN09780 
IE  #RSnLUTE  VALUE  OF  F(Z)  .LE.  ePS  (INPUT)  ZAN09790 
2ND  STOPPING  CRITERION.  A BOOT  IS  ACCEPTED  ZAN09800 
IF  TWO  SUCCESSIVE  APPROXIMATIONS  TO  A GIVEN  ZAN09B10 
ROnT  AGREE  IN  THE  FIRST  NSIG  DIGITS.  (I NPUT) Z AN09820 
NOTE.  |F  E TTHFR  CR  BOTH  O F thf  STOPPING  ZAN09830 

OB  I TF° I A A»F  FULFILLED,  THF  RCCT-  IS  ZAN09840 

AOCEOTED.  ZAN0985O 

THF  NUMREP  CF  KNOWN  BOCTS  WHICH  MUST  RF  STOR EDZ AN09860 
T M X(  U,...,X(KN),  PRIOR  TO  ENTFY  to  ZANLYT  ZAN09870 
THE  NUMPFP  of  initial  guesses  provided.  THESE 
GUESSES  MUST  RF  STORED  in  X(KN+1),..., 

X(  KN-fNGHESS)  AND  NGUESS  MUST  BE  SET  EQUAL 
TO  ZFP.O  IF  Nn  GUESS  FS  AP  F PROVIDED.  (INPUT) 

THF  NUMREP  OF  NEW  ROOTS  TO  RE  FOUND  BY 
ZANLYT  (INPUT) 

LONG-WOPD  COMPLEX  VECTOB  ARRAY  OF  LENGTH 
.GE.  3*(KN*N).  X(1 ) ,.. .,X(KN)  ON  TNPUT 
MijST  CONTAIN  ANY  KNOWN  POOTS.  X(  KN*  1 ) , . . . , 

X(KN*N)  ON  INPUT  MAY,  AT  THE  USFP'S  OPTION, 

CONTAIN  INITIAL  GUESSES  FOR  THE  N NEW 
ponr$  WHICH  APE  TO  BE  COMPU^FD.  ON  OUTPUT, 

XIKM*!  ),...,  X(K.N*N)  CONTAIN  E!  THFb  A P"OT 
CORRpry  TO  WITHIN  A CONVERGENCE  CPITFPPN 
OP  THE  V ALU E(12 345678. 123456780+0, 12345678. 
123*56780*0)  INDICATIVE  OF  A FAILURE  T0 
ACHIFVc  THE  SPFCIFItD  CONVERGENCE  FOR  THAT 
PPOT.,  «4Y  X ( K N ♦ J ) . IN  THE  LATTEF  CASF,  tH£ 

MOST  RFC ENT  APPROXIMATION  TO  XIKN+J)  IS 
AV  AM  ABLE  IN  XIISUB),  WHFRF  I SU«  = 2*  ( KN«-N  )+J 
THF  MAXIMUM  ALLOWABLE  NUNPFR  OF  iterations 
PFR  prnT  (INPUT) 

AN  INTEGFP  VFCTTP  CF  LENGTH  .GE.  KN+N.  ON 
OUTPUT  TNFFR(J)  contains  the  number  OR 
TT'PP  6T I O')  S USED  IN  FINDING  THE  J-TH  ROOT 
WHEN  CONVERGENCE  WAR  ACHIEVED.  IP 
CONVERGENCE  WAS  NOT  ORTATNFD  IN  ITMAX 
ITccatiovs,  IMFFR(J)  WILL  CONTAIN  ITMAXU 


C 

o. 

r 


IEP 


PRECISION 

PEQ'D  I M$L  POUT  INFS 
AUTHOP/IMPLEMENTOP 
L ANGUAGF 


(OUTPUT) 

FPPno  parameter  (OUTPUT) 

WARNING  RRROP  = 32  + N 

M * 7 FAitUBF  tp  CONVERGE 
ITPPATIONS  P3P  rwF  OF  THF 
BE  ROUNO 
DOUBLE 
UPBtst 

o.  G.  JOHNSCN/L.  L.  WILLIAMS 

FORtpan 


WITHIN  ITMAX 
(M)  NFW  BOOTS  TC 


latest  bfvision 


SEPTEMBER  1,  1971 

65 


ZAN09880 
Z AN0989D 
ZAN0990D 
Z AN09910 
ZAN09920 
ZAN09930 
Z AN09940 
Z AN09950 
Z AN09960 
ZAN09970 
ZAN0998D 
ZAN09990 
ZANIOOOO 
ZANI0010 
ZAN10020 
ZAN10Q30 
Z * Ml 0040 
ZAN10050 
Z4NI0060 
ZAN1O07O 
ZAN10083 
ZAN1039D 
ZAN10100 
ZAN10110 
ZAN1012O 
Z AN  101 30 
Z AN  1 0 1 40 
ZAN1C150 
Z AN  101  60 
ZAN101  70 
ZAN1018" 
Z AN 101 90 
ZAN10200 
Z AN  102 10 
ZAN1Q220 
ZAN10230 
Z AN10P40 
Z AM  02  50 
, Z A N1  0?  6” 
Z AN  10? 70 


r.  K • -yrli.'Xl-U/  -M*  **-Z-*~ 


SUBROUTINE  ZANLYT 
C0MPLFX*16 


DOUBLE  PRECISION  OZ,Fi»?,FP 

OI«ENSION  INFFP(l) 

IEP  « 0 

ONE  - (1.00*00,0.00*00) 

FPS1  ■ 10.00*00**(-NSTG) 

TCONJ  « 0 
IROMB  « 0 

( 

MR  1 - KN*1 
MR 2 ■ KN*N 
1ST  APT  ■ MB 2*1 
MP3  « MB1*NGUESS 
DO  ? I * MPG.MB2 

2 X ( 1 1 « (0.0*0,0«D*0) 

L * MBl 

IF  (KN  .EQ.  0)  GO  TO  5 
DO  3 T * 1 » KN 
INFER(I)  « 0 
TTFMP  a MR2+I 
XHTEMD)  * XCI  ) 

ITEMP  a MR  2*ITEMP 

3 X(TTEMP)  * X(I) 

5 JK  a 0 

oz  - CDABsmm 

IF  (OZ  .LE.  l.OD-15)  GO  Tn  25 

f 

10  PT  a ( .90*00, 0.00*00)*X(L) 
ASSIGN  15  TO  NM 
GO  T0  135 
15  XO  a FPPT 

PT  a u.  10*00, 0.00+00)*X(L) 
ASSIGN  20  TO  NN 
GO  TO  135 
20  XI  = FPPT 
H a X(L)-RT 
PT  a X ( L ) 

ASSIGN  40  to  NN 
GO  Ti  135 


<F,FPS,NSIG,KN,NGUESS,N,X,  ITMAX, INFER, IEP ) 
XI 1)  ,f'NE,D,DD,DFN,D!,FPRT,FPT, 
H,RT,T1,T2,T3,TEM,X0,X1,X2,BI ,F,XX 
0Z,FPS,FPS1 
INFFP ( 1) 


SET  NUMBER  OF  ITERATIONS 


POCT  FSTJMATC  NO*  equal  TO  ZERO 


ROOT  ESTIMATE  EQUAL  T0  ZFPC 


25  RT  = -ONE 

ASSIGN  30  TO  NN 
GC  Tn  135 
30  XO  = FPRT 
»T  a ONE 
ASST3N  35  TO  NN 
GO  TO  1*5 
35  XI  = FPRT 

P’  a (0.00*00,0.00*00) 
H * -ONE 
ASSTGN  40  TP  NN 
3P  to  135 
40  X2  = PPPT 

45  D = (-0.50*00,0.00*00) 

50  00  = ONF  ♦ 0 
▼1  a X0*0*0 
T2  a X1*00*DO 


RFGIN  main  ALGORITHM 


ZAN10280 
Z AN10290 
ZAN10300 
ZAN1031P 
Z AN10320 
Z AN10330 
ZAN10340 
ZAN10350 
ZAN10360 
Z AN10370 
ZAN103B0 
Z AN10390 
Z AN10400 
Z AN) 0410 
ZAN10420 
Z AN10430 
Z AN10440 
ZAN10450 
Z AN10460 
Z AN10470 
ZAN10480 
Z AN10490 
Z &N1O500 
Z AN10510 
ZAN10520 
ZAN1053O 
Z AN  10540 
ZAN10550 
ZAN10560 
Z AN10570 
ZAN10580 
Z AN10590 
ZAN10600 
Z A Ml  06 10 
Z 4N10620 
Z AN! 0630 
ZANIC640 
Z AN  10650 
ZAN10660 
ZAN10670 
Z AN10680 
Z AN10590 
ZAN10700 
Z AN10710 
ZAN10720 
ZAN10730 
ZAN10740 
Z AN10750 
Z AN10760 
Z AN10770 
ZAN10780 
Z AN10790 
ZAN1080D 
ZAN10810 
Z AN10820 
ZAN10830 
ZAN1P840 
Z AN10850 
l AN10860 

zani'-rzo 

Z A N1  08  80 


55 

60 

65 

70 

75 


r 

r 


80 

85 

00 

05 


XX  « X2*D0 

-3  = X2*D 

«!  « T1-t2+XX*T3 

OEM  * RI*R!-(4.00*-0,o.0«-P)  *(  XX*Ti-T  3*(T2-XX)  ) 

»JSF  OENCMNATPR  r F MAXIMUM  4MPI  ITUOt 

ri  * roso«T(riFM) 

72  * «I  ♦ T1 
n > - t| 

OZ  * rOARS(T2)  - rnAPSC'3) 

IF  (QZ  .OF.  m r,n  tP  60 
Oc M * T3 
r,r  T0  65 
OFN  * T? 

Tp  ST  FfP  ZF-H  OFMOMtKATOP 

QZ  = rOAP5(OEN) 

IF  HZ  .OT.  1.0-15)  on  -r  75 

pfM  s ns|F 

or  * ((-?. "pop, o.oruoo>*xx)/oeN 

h * or  * w 

PT  = RT  ♦ W 

OMFCK  CON'V  FFGFMf  E Q F THF  FIRST  KINO 

OZ  = 'OABSIH/PT) 

IF  (OZ  .L  E.  FPSI)  O0  rn  l0p 
AFC TpN  85  rr  mm 
op  tp  135 

r>Z  * C DAPS  ( FPRT  l-CP APS  (X2* ( 1°  .3  00*0  .O0r  ) ) 

IF  (3Z  .LT.  0.P*0)  On  tp  05 

TA  KF  OEMFppL  Af’irf!  T?  TNOIJOF 
OCNVFPOFMOC 

DT  = 0T*(  0.50*00,0.0^+00) 

H = H*(Q .50*00, 0.00+00) 

PT  = BT_8 
OP  Tp  135 

XO  = XI 
XI  * X2 
X?  a FPPT 
0 * 01 
",n  TO  50 


A FPPT  HAS  PFFN  FPIINC 


100 

105 


FRT  = F(P'r) 

XI  L ) = RT 

TTCMp  a MO?*l-lBOMP 
XUTFMP)  a BT 
1 r p mo  - M32+MP7+L 
X ( !TFMP ) = ot 


,0 


phfck  tp  ccf  tf  rpMPLFx-rrr’juoxTF 
TS  USC  A Bn0^ 

10.r*O*PDABS( FPT) ) GO  tf  115 


\ 


110 


115 

120 

125 


TF  (CD  APS  ( F(  D<TMJ0(  X (L  ) ) ) ) 

OZ  = CD»RMX<U-  orPMjrixd  ) ) ) 

IF  ( TC°NJ  .ME.  0 .PP.  OZ  . L T . 1.00-B)  GO  rr  ))5 
tsta»t  = 1*2 
T V S F p T = L + l 

on  un  INSFPT  a TS'rACT,M32 
X(  INSERT)  = X(INSC°1I 
T V)  c F P 1 a I N S F B ',‘ 

X ( L +1  ) = OrnNjr,(x(i  )) 

T-OMJ  = l 
00  T-1  120 
irpvj  = ry 
:,nMTT  Mljf 
TNFF3((_)  a JK 


?, AN1C89P 
Z A'll  0900 

zAnio?n 

Z '.Ml  0920 
Z4N109OP 
ZAN10940 
Z AN 109 50 
ZANIP96D 
ZAM  p 9 7 n 
Z AN1098P 
ZAN1099D 
ZAN1  JPP.? 
Z AN1 1010 
ZAN11020 
Z AN  110  30 
Z ANT  )040 
ZAM105D 
Z AN  11060 
Z AM  10  70 
ZAN11080 
ZAN11090 
ZAN11100 
ZAN11110 
ZAN11120 
ZANH130 
ZAM  11  40 
Z AN  111  50 
Z AM  1160 
Z AN1 T 1 70 
Z ANtllBO 
Z4N1119D 
Z ?N1 1200 
ZANU71P 
Z AN1 1220 
Z AN  11230 
ZAN11240 
ZAN1125P 
ZANU260 
ZAM11270 

Z AN1 1280 
Z AM  1 2 90 
Z AM  1303 
Z AM  1 31  o 
ZANU  320 
Z AM  i 330 
ZAN1 1 340 
7 AN) ) o 50 
ZAM  1360 
ZAN11370 
Z AMI  13  80 
7.ANU390 
ZAM  1 400 
Z AMI  14  ID 
ZAN1142n 
Z AMT  1430 
ZAN11 44D 
Z AM  14  50 
Z ANl T 460 

ZAM.  I 4 70 
Z a N 1 1480 
ZAM149P 
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•a?!"1"!'  "J 


Mi"  i . 

Z_3 


IlMIl,  V ..M  M |l  fill)) 


L « 1*1 

1 1 * J 1 1,  5"0 

IF  (L  .IF.  M97)  G"  ”0  5 

Z * M 1 5 1 0 

c 

’FTUF'f  Tn  C A l L J NC  PProoftM 

Z4MU52G 

130 

jo  TO  185 

£ * M 1 15  30 

135 

JK  * JKM 

Z 8 Nil  54', 

T F (JK  .GT.  TTM1XI  r,n  Tn  \ 90 

Z4HH550 

140 

FPT  = F ( P T j 

Z *N! 1560 

PDPT  = FPT 

Z 4 n*  1 1570 

r 

rrc-  re  rrc  jp  r'ppT 

i <;  BEING 

Z A N 1 1 5 rtr 

r 

Oftfp  HI \'CD 

Z 8\J11  590 

IF  (L  .FO.  n 00  TO  160 

Z ft  Ml  15  00 

if  (l  #lf . iromr+ij  r-n  t"  i6^ 

Z » Nil  610 

r 

roMPij’-P  OFN^MIMftTfp  FOR  ^PCIFIFD 

Z ft\'l  1 620 

f 

F'JNCTIr'N 

Z ft m 1 1 630 

145 

LT  HUP  = MR2  *l-T  ROmr-i 

7.  AN  1 1640 

00  150  I = LST4R‘r,lTvllP 

Z4M165C 

TCM  X,  RT  - X(T  ) 

Z i'll  1 660 

OZ  a r08P5(TFMJ 

Z*N1 1670 

TF  (OZ  . L T « 5.00-1  5»  O'*-  ‘rn  175 

Zft‘JU480 

150 

FpOT  * FPPT/TFM 

Z ft  MIJ  6 9,o 

FHFFK  CFfiVFrOFNFF  OF 

THC  S FC  “’NO  K INI) 

Z IN  11  ZOO 

160 

OZ  = C0ABS( F3T| 

Z ft.  M 1 1713 

IF  (OZ  .OF.  FPC)  G0  Tn  IT? 

ZftMU720 

165 

oz  * :o4rs(fpdtj 

1C  Nil  7 30 

if  (OZ  .IT.  EPF  ) or  To  105 

Z4N11740 

170 

Or  Tr  NN,(15, 70,30, 35,40,86) 

Z ftMl 1750 

175 

PT  « ?T  ♦ ( 1 .OOOOOimo.O.GfHO) 

Z 'Ml  1 760 

G"1  TO  135 

Z ft  Ml  1 J 7? 

r 

WftFNING  FP  P r 3 , ITMAX 

r MA  X T M 1 JM 

Z INI  1780 

180 

TFO  a 37 

Z 8N 1 1 790 

INFFO(L)  = T t Mf  x *■  1 

Z ft  Ml  1 

tr.omr  = jrohr  + 1 

z ini  mo 

X ( L ) = (17345678. 1 ?345678n  + n, 1 ? 3456 78.  1 23456780*0) 

Z ft  M 11  -1 2 0 

JTFMP  = MR2  «■  MR2  + 1. 

Z ft  N 1 1.830 

X( I TFMP ) = PT 

ZftNl)84n 

L = L ♦ ! 

Z ftMl 1850 

IF  (L  .IE.  mrp  ) GO  Tr'  5 

Z ft  N 1 1 ft  60 

185 

IF  ( IF*  .FO . 0)  or  7n  9005 

Z ft N 118  70 

9000 

OOMT  IMIJF 

Z ftMl  l flp'- 

:*ll  UEPTST (TF0,*  ZCWLYT’  ) 

Z 'Nil  690 

9005 

3 = T(IPN 

Z ft'1 1 1900 

PMO 

Z ftMl  1 91  o 
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1 o 


• T'R.V'UWP'H  JBIM.'I 


1 


(JFf>l  P?? 

o 'J  F ® 1 1 9 3 ^ 

r FUNr  TT  9N  - eopre  w£s«;f(;F  GFNEF'TTPN  ijcr  \1  94C 

r ij$4gp  - tai  l hcct«;t(  jF,  t mEXX,  j ijfru950 

r PARSMfTpqs  IF®  - FPfpp  o^ofMprfF.  tvpf  * »•  WHFP  F ljec'M96r 

r TyoFa  i?g  imoltFc  tpp«!'i7I.  frohc  jc&\io7r, 

r 64  JMPL’EF  WANING  WITH  fjX  IJFP119P0 

r 3?  Tmpltfc  Wf  c N Tf’G  ijpo  1. 1 4 99 

N = cpcrc  CTDF  PFl  fv.rkt  T"  CM  LING  F pjt  j vf  j cc i 0903 
N4MFXX  - M*MC  IF  TWF  f £ l l TNG  on,JMNp  JEP1221C 

r AUTHFT/T  MPLFMENTfh  _ of  nFe  8VCMPSFN  JFP 12020 

r LANGUAGE  - fortps^  JE°12030 

r jff 1294c 

C L»TF?T  PFVIG!?^  - JMJII»oy  19,  1971  UCP1?:050 

r GF"  12960 

StJFPHjTlNF  UFBTPT(tfp,mcmcj  tjco  1 2270 

r J Fa \ 20  80 


tUMENS  THN 

TTYo  ( 5,  4)  , 1 P JT{  4 ) 

UF  R l 209r> 

TSjTpr,FR*2 

M4MC (H  ) 

'JPP  l?100 

TNTFGCP 

WfcON,W«OF,TFPM,pQTNiTp 

JE°1 21  1C 

cOJT\MLFN 

CE 

(»PTT(1  |,W«M,(  I R IT  ( 2),  W1FF  »,  ( IRI7  (3)  UFF*) 

11 F f 1 2 1 22 

n»Tf 

t TYP 

/'W't'*"  ,*  TMG  • , • • , • • , 1 • , 

IIEP  12’  30 

* 

•W7PM* , 'TNO ( • , • WITH* , • MX',')  ', 

UEP l 2 1 40 

* 

• TFom<  , • TMU'»  • • , • ' • , 

'Jpo  1 2 I 5^ 

* 

•NPN-* ,'nFF I* , • NFO  •/, 

JFP  1 21  60 

* 

T®  JT 

/ 32.64,128,0/ 

'JPF12179 

r>*T* 

PPJNTp  / 6/ 

IJFD  1 ?l«9 

1co?sTfd 

JE0?  21 9C 

IF  Uco? 

.GF. 

WARN)  Gn  Tr  5 

IJFC  1 22  nr' 

r 

< 

N9K-0FPTMP0 

1JPF  12210 

IFO 1 =4 

'JPR  1 2220 

Go  tp  ?f> 

•JEP  1 2239 

5 

IF  ( T cp  ? 

.LT. 

tfpm)  r.n  Tn  m 

yfc i ??40 

r 

-rF  p HT  \A  l 

'JpR  1 22  50 

T CP 1 = 3 

ijf’  l ??fer 

3°  T9  29 

UCR  J 2?  70 

10 

TP  (IFF? 

• IT. 

WAFF)  Gn  -n  15 

1 IPR  l ?£  80 

r 

(WITH  FIX) 

ijco  i2?9o 

!P5  l=? 

ifoi ?Tnr 

ro  T7  20 

UE°  12  310 

r 

w A PH  FNG 

Jpo 13?  20 

15 

T P F 1 * 1 

JFS1 2? 30 

r 

c x Tp  ft  r t I..II 

1.1  FA  1 ,-M40 

29 

T F°  2=1 CF  2 

-I  PI  T | 1 FP 1 ) 

•jro  1 23  50 

r 

or  j fiT  fdppc  vp<;c(  r,c 
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